Skip to main contentIBM Quantum Documentation Preview
This is a preview build of IBM Quantum® documentation. Refer to quantum.cloud.ibm.com/docs for the official documentation.

Long-range entanglement with dynamic circuits

Usage estimate: 90 minutes on IBM Sherbrooke (NOTE: This is an estimate only. Your runtime may vary.)

Background

Long-range entanglement between distant qubits on a quantum processor can be a challenging task for devices with limited qubit connectivity. This tutorial shows three different ways that can be used to generate long-range entanglement between qubits on a line, at varying distances between each other:

  • a unitary-based implementation, which uses SWAP operations to reduce the distance between distant qubits and entangle them directly,
  • a measurement-based implementation with post-processing, which discards some amount of information to generate the desired entangled state, and
  • a measurement-based implementation with dynamic circuits, which uses measurement and feedforward of information during the quantum computation to entangle the qubits. In particular, the results show the effectiveness of dynamic circuits in generating long-range entanglement between two unconnected qubits at utility scales. The notebook uses the ideas and results from [1] by Elisa Bäumer et al.

Requirements

Before starting this tutorial, ensure that you have the following installed:

  • Qiskit SDK 2.0 or later, with visualization support ( pip install 'qiskit[visualization]' )
  • Qiskit Runtime ( pip install qiskit-ibm-runtime ) 0.37 or later

Setup

import random
from typing import List, Dict, Union, Optional, Callable
 
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import clear_output
 
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit import Gate
from qiskit.circuit.classical import expr
from qiskit.circuit.library import XGate
from qiskit.primitives import BitArray
from qiskit.providers.backend import BackendV2 as Backend
from qiskit.quantum_info import Pauli, PauliList
from qiskit.result import marginal_counts
from qiskit.transpiler import CouplingMap, InstructionDurations
from qiskit.transpiler.passmanager import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import Batch, SamplerV2 as Sampler
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime.transpiler.passes.scheduling import (
    DynamicCircuitInstructionDurations,
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)

Step 1: Map classical inputs to a quantum problem

In this tutorial you will run a gate teleportation circuit in three different setups, where you always assume a line of n qubits (for varying n with n-2 empty ancillas in the middle and a CNOT to apply between the two ends):

  • Unitary-based implementation swapping the qubits to the middle
  • Measurement-based implementation with post-processing
  • Measurement-based implementation with dynamic circuits

For each implementation you measure the average gate fidelity to compare among the different implementations. For details on how the average gate fidelity is calculated, see the Appendix.

Experimental setup

The experiments in this notebook use a 1-D line of qubits with a coupling map that ensures that no shortcuts can be taken.

Define 1-D line

First, set up a line of qubits through the machine that you intend to use such that you avoid broken qubits or areas with high readout errors. To do this, examine the calibration data (which can be found online or via the command plot_error_map(backend)). For example, suppose that you use ibm_sherbrooke and that you need to avoid, for example, qubits 20 and 71 as well as qubits 56 and 73. One such qubit line would be:

A coupling map of ibm_sherbrooke showing a subset of qubits connected by a single continuous line

We will describe the line as a simple list of integer indices. In this tutorial, we will choose a random qubit line using the following function.

def get_1d_qubit_line(
    coupling_map: CouplingMap, start_qubit: int = 0
) -> List[int]:
    """
    Use random search to find a 1D line, repeating until we find a line of sufficient length.
    This is very inefficient, but it's fine for this demo.
    """
 
    def get_random_line(
        coupling_map: CouplingMap, start_qubit: int
    ) -> List[int]:
        """
        Do a random, self-avoiding walk on the coupling map to get a 1D line.
        """
        edge_list = list(coupling_map.get_edges())
        path = [start_qubit]
        while True:
            # Get edges connected to the current qubit that don't connect to a qubit we've already visited
            current_qubit = path[-1]
            previously_visited_qubits = path[:-1]
            candidate_edges = [
                edge
                for edge in edge_list
                if current_qubit in edge
                and not any(
                    qubit in previously_visited_qubits for qubit in edge
                )
            ]
            if candidate_edges == []:
                return path
            next_edge = random.choice(candidate_edges)
            edge_list.remove(next_edge)
            next_qubit = next(
                qubit for qubit in next_edge if qubit != path[-1]
            )
            path.append(next_qubit)
 
    # Now repeat the random walk many times to find a long line
    return max(
        [get_random_line(coupling_map, start_qubit) for _ in range(1000)],
        key=lambda x: len(x),
    )

The following cell shows an example of such a line.

from qiskit_ibm_runtime.fake_provider import FakeSherbrooke
 
qubit_line = get_1d_qubit_line(FakeSherbrooke().coupling_map)
print(f"Found line of length {len(qubit_line)}:")
print(" → ".join(str(q) for q in qubit_line))

Output:

Found line of length 76:
0 → 14 → 18 → 19 → 20 → 33 → 39 → 38 → 37 → 52 → 56 → 57 → 58 → 71 → 77 → 78 → 79 → 91 → 98 → 99 → 100 → 101 → 102 → 92 → 83 → 84 → 85 → 73 → 66 → 65 → 64 → 54 → 45 → 44 → 43 → 34 → 24 → 25 → 26 → 16 → 8 → 9 → 10 → 11 → 12 → 17 → 30 → 31 → 32 → 36 → 51 → 50 → 49 → 55 → 68 → 69 → 70 → 74 → 89 → 88 → 87 → 93 → 106 → 107 → 108 → 112 → 126 → 125 → 124 → 123 → 122 → 121 → 120 → 119 → 118 → 110

Set primary parameters

In this section are definitions for some common parameters that you will use later. You'll need to specify these parameters for a particular backend. In order to do so, you will need an account on IBM Quantum™ Platform. More details on how to initialize your account can be found in the documentation.

def coupling_map_from_qubit_line(
    coupling_map: CouplingMap, qubit_line: List[List[int]]
) -> CouplingMap:
    """
    Modify the full coupling map to force linearity in the qubit layout
    """
    new_coupling_map = []
    # Iterate the line pair-wise and append edges that contain both qubits
    for qubit, next_qubit in zip(qubit_line, qubit_line[1:]):
        edge = next(
            edge
            for edge in coupling_map.get_edges()
            if qubit in edge and next_qubit in edge
        )
        new_coupling_map.append(edge)
    return CouplingMap(new_coupling_map)
 
 
# Set up access to IBM Quantum devices
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=127
)
 
# Set qubit line and coupling map
QUBIT_LINE = get_1d_qubit_line(backend.coupling_map)
COUPLING_MAP_1D = coupling_map_from_qubit_line(
    backend.coupling_map, QUBIT_LINE
)
MAX_POSSIBLE_QUBITS_BTW_CNOT = len(QUBIT_LINE) - 2
 
# Use this duration class to get appropriate durations for dynamic
# circuit backend scheduling
DURATIONS = DynamicCircuitInstructionDurations.from_backend(backend)
 
print(f"Backend is: {backend.name}")
print(f"Found line of length {len(QUBIT_LINE)}.")
print("Using coupling map:\n", COUPLING_MAP_1D)
print(
    f"Maximum number of qubits between CNOT for {backend.name} is {MAX_POSSIBLE_QUBITS_BTW_CNOT} with the given qubit line."
)

Output:

Backend is: ibm_brisbane
Found line of length 76.
Using coupling map:
 [[14, 0], [14, 18], [18, 19], [20, 19], [21, 20], [21, 22], [15, 22], [4, 15], [4, 5], [6, 5], [6, 7], [7, 8], [16, 8], [16, 26], [27, 26], [28, 27], [28, 29], [30, 29], [30, 31], [31, 32], [32, 36], [36, 51], [50, 51], [50, 49], [48, 49], [48, 47], [46, 47], [46, 45], [44, 45], [43, 44], [42, 43], [42, 41], [41, 53], [53, 60], [60, 61], [62, 61], [62, 72], [81, 72], [81, 82], [82, 83], [84, 83], [85, 84], [85, 86], [86, 87], [93, 87], [93, 106], [105, 106], [105, 104], [111, 104], [122, 111], [122, 121], [121, 120], [120, 119], [118, 119], [117, 118], [116, 117], [116, 115], [114, 115], [114, 109], [109, 96], [97, 96], [97, 98], [98, 91], [91, 79], [79, 78], [77, 78], [77, 71], [58, 71], [57, 58], [56, 57], [52, 56], [52, 37], [37, 38], [39, 38], [40, 39]]
Maximum number of qubits between CNOT for ibm_brisbane is 74 with the given qubit line.

Next, set the global scope of the experiment. These variables can be used in each circuit type or can be set individually in each experiment that will override these globals.

# Set which Paulis to sample (default is all 16 combinations that have a non-zero expectation)
SAMPLES = list(range(16))
 
# Level of optimizations the transpiler uses: There are 4 optimization levels ranging from 0 to 3,
# where higher optimization levels take more time and computational effort but may yield a more
# optimal circuit. level 0 does no explicit optimization, it will just try to make a circuit
# runnable by transforming it to match a topology and basis gate set, if necessary.
# Level 1, 2 and 3 do light, medium and heavy optimization, using a combination of passes, and by
# configuring the passes to search for better solutions.
OPTIMIZATION_LEVEL = 1
 
# Set to True to use dynamical decoupling
USE_DYNAMIC_DECOUPLING = True
 
# Default dynamic decoupling sequence if dynamic decoupling is used
DD_SEQUENCE = [XGate(), XGate()]
 
# Set the number of repetitions of each circuit, for sampling.
# The number of qubits between the control and target are grouped into blocks
# of length 4. The provided min and max number of qubits will be modified to
# align with these block sizes.
SHOTS = 1000
 
# The min number of qubits between the control and target qubits on line
MIN_NUMBER_QUBITS = 0
 
# The max number of qubits between the control and target qubits on line
# The max for MIN_NUMBER_QUBITS is len(QUBIT_LINE) - 2
# max_number_qubits must satisfy MAX_NUMBER_QUBITS - MIN_NUMBER_QUBITS = 3 (mod 4)
# at this point. This is just to make things easier to break jobs up. Not a real limitation.
MAX_NUMBER_QUBITS = 63
assert (
    (MAX_NUMBER_QUBITS - MIN_NUMBER_QUBITS) % 4 == 3
), "MAX_NUMBER_QUBITS must satisfy MAX_NUMBER_QUBITS - MIN_NUMBER_QUBITS = 3 (mod 4)"
assert (MAX_NUMBER_QUBITS + 2) <= len(
    QUBIT_LINE
), "MAX_NUMBER_QUBITS must satisfy MAX_NUMBER_QUBITS + 2 <= len(QUBIT_LINE)"

Monte Carlo state certification

Next, you utilize the methods prep_P_ij_conj and meas_P_kl to prepare the circuits for Monte Carlo state certification. The method can be used to calculate the average state fidelity of the state prepared by a circuit without having to estimate the full quantum state (for example with state tomography or similar technqiques). For more details see the Appendix. Monte Carlo state certification requires us to prepare the complex conjugate of random product of eigenstates of local Pauli operators, corresponding to the Paulis Pi∗P_i^* and Pj∗P_j^* - and then to measure the the final state in different Pauli bases, corresponding to the Pauli operators PkP_k and PlP_l. You can do this with the following two methods: the prep_P_ij_conj method prepares a list of circuits so that the control and target qubits are in the eigenstates of Pi∗P_i^* and Pj∗P_j^*, respectively; and the meas_P_kl method prepares a list of circuits so that the control and target qubits are measured in the PkP_k and PlP_l bases, respectively.

def prep_P_ij_conj(
    circuits: List[QuantumCircuit], P_prep: PauliList
) -> List[QuantumCircuit]:
    """Prepare circuits with the possible complex conjugates of the eigenstates of given Pauli operators
 
    The first and last qubits are prepared in one of the four possible eigenstates of P_i^* and P_j^*
    respectively. The resulting collection of circuits covers all four possibilities.
 
    Assumes that circuits have qubits 0, ... ,n+1 where the long range NOT is between qubit 0 (control) to n+1 (target)
 
    Arg:
        circuits (List[QuantumCircuit]: List of 4 Quantum Circuits with at least two data qubits
        P_prep (PauliList): A pair of single qubit Paulis with the first Pauli referred to as P_i and the second as P_j
 
    Returns:
        List[QuantumCircuits]
    """
 
    pauliX = Pauli("X")
    pauliY = Pauli("Y")
    for p in range(2):
        for q in range(2):
            qc = circuits[2 * p + q]
            if p == 1:
                qc.x(0)
            if q == 1:
                qc.x(-1)
            for i in range(2):
                if P_prep[i] == pauliX:
                    qc.h(-i)
                if P_prep[i] == pauliY:
                    qc.h(-i)  # Change basis to initialise in Y^* basis
                    qc.s(-i)  # i.e. Apply SH = (HS^\dagger)^*
            circuits[2 * p + q] = qc
    return circuits
 
 
def meas_P_kl(
    circuits: List[QuantumCircuit], P_meas: PauliList
) -> List[QuantumCircuit]:
    """Prepare circuits so that the final state can be measured in different Pauli bases with a Z measurement
 
    The first and last qubits are the qubits that the given operator P_meas will operate
 
    Arg:
        circuits (List[QuantumCircuit]: List of 4 Quantum Circuits with at least two data qubits
        P_meas (PauliList): A pair of single qubit Paulis with the first Pauli referred to as P_k and the second as P_l
    """
    pauliX = Pauli("X")
    pauliY = Pauli("Y")
    for p in range(2):
        for q in range(2):
            qc = circuits[2 * p + q]
            for i in range(2):
                if P_meas[i] == pauliX:
                    qc.h(
                        -i
                    )  # Change of basis to measure X by measuring in Z basis: i.e., apply H
                if P_meas[i] == pauliY:
                    qc.sdg(
                        -i
                    )  # Change of basis to measure Y by measuring in Z basis
                    qc.h(-i)  # i.e. apply HS^dagger
            qc.measure(0, 0)
            qc.measure(-1, 1)
            circuits[2 * p + q] = qc
    return circuits

Unitary-based implementation swapping the qubits to the middle

First, examine the case where a long-range CNOT gate is implemented using nearest-neighbor connections and unitary gates. In the following figure, on the left is a circuit for a long-range CNOT gate spanning a 1D chain of n-qubits subject to nearest-neighbor connections only. On the right is an equivalent unitary decomposition implementable with local CNOT gates, circuit depth O(n).

Two circuit diagrams, one showing a CNOT acting on non-neighbor qubits, and the other showing a sequence of neighbor-only CNOTs that achieve the same behaviour.

The circuit on the right can be implemented as follows:

def CNOT_unitary(
    qc: QuantumCircuit, control_qubit: int, target_qubit: int
) -> QuantumCircuit:
    """Generate a CNOT gate between data qubit control_qubit and data qubit target_qubit using local CNOTs
 
    Assumes that the long-range CNOT gate will be spanning a 1D chain of n-qubits subject to nearest-neighbor
    connections only with the chain starting at the control qubit and finishing at the target qubit.
 
    Assumes that control_qubit < target_qubit (as integers) and that the provided circuit qc has |0> set
    qubits control_qubit+1, ..., target_qubit-1
 
    Args:
        qc (QuantumCicruit) : A Quantum Circuit to add the long range localized unitary CNOT
        control_qubit (int) : The qubit used as the control.
        target_qubi (int) : The qubit targeted by the gate.
 
    Example:
 
        qc = QuantumCircuit(8,2)
        qc = CNOT_unitary(qc, 0, 7)
 
    Returns:
        QuantumCircuit
 
    """
    assert target_qubit > control_qubit
    n = target_qubit - control_qubit - 1
    k = int(n / 2)
    qc.barrier()
    for i in range(control_qubit, control_qubit + k):
        qc.cx(i, i + 1)
        qc.cx(i + 1, i)
        qc.cx(-i - 1, -i - 2)
        qc.cx(-i - 2, -i - 1)
    if n % 2 == 1:
        qc.cx(k + 2, k + 1)
        qc.cx(k + 1, k + 2)
    qc.barrier()
    qc.cx(k, k + 1)
    for i in range(control_qubit, control_qubit + k):
        qc.cx(k - i, k - 1 - i)
        qc.cx(k - 1 - i, k - i)
        qc.cx(k + i + 1, k + i + 2)
        qc.cx(k + i + 2, k + i + 1)
    if n % 2 == 1:
        qc.cx(-2, -1)
        qc.cx(-1, -2)
    qc.barrier()
    return qc

The build_circuits_uni method therefore builds a list of circuits to run with different Paulis Pi,Pj,PkP_i, P_j, P_k and PlP_l in order to estimate the state fidelity with the Monte Carlo state certification method.

def build_circuits_uni(n: int, samples: List[int]) -> List[QuantumCircuit]:
    """Builds the unitary circuits needed to estimate the average gate fidelity
 
    Args:
        samples (List[int]): Which of the 16 Paulis with non-zero expectation value to prepare and measure
        n (int): Number of qubits between the control and target of the CNOT
    """
    circuits_all = []
    # 16 Paulis with non-zero expectation value to prepare and measure
    P_lkji = PauliList(
        [
            "IIII",
            "XIXI",
            "IZIZ",
            "XZXZ",
            "YZYI",
            "ZZZI",
            "YIYZ",
            "ZIZZ",
            "XXIX",
            "IXXX",
            "XYIY",
            "IYXY",
            "ZYYX",
            "YYZX",
            "ZXYY",
            "YXZY",
        ]
    )
    for sample in samples:
        P_prep = P_lkji[sample][0:2]
        P_meas = P_lkji[sample][2:4]
        circuits = [QuantumCircuit(n + 2, 2) for i in range(4)]
        circuits = prep_P_ij_conj(
            circuits, P_prep
        )  # Prepare circuits in eigenstates P_i^* and P_j^*
        circuits = [
            CNOT_unitary(circuit, 0, n + 1) for circuit in circuits
        ]  # Add long range CNOT
        circuits = meas_P_kl(
            circuits, P_meas
        )  # Prepare circuits to measure the final state in P_k and P_l bases
        circuits_all += circuits
    return circuits_all

For example:

# Sample circuit
n = 6
sample = [11]
test_circuits = build_circuits_uni(n, sample)
test_circuits[3].draw("mpl", fold=-1)

Output:

Output of the previous code cell

Measurement-based implementation with post-processing

Next, examine the case where a long-range CNOT gate is implemented using nearest-neighbor connections of a measurement-based CNOT with post-processing. In the following figure, below on the left is a circuit for a long-range CNOT gate spanning a 1D chain of n-qubits subject to nearest-neighbor connections only. On the right is an equivalent implementable with local CNOT gates with measurements, and which requires post-processing.

Two circuit diagrams showing a single CNOT operation before and after post processing. Both circuits are equivalent.

The circuit on the right can be implemented as follows:

def CNOT_postproc(
    qc: QuantumCircuit,
    control_qubit: int,
    target_qubit: int,
    c1: Optional[ClassicalRegister] = None,
    c2: Optional[ClassicalRegister] = None,
    add_barriers: Optional[bool] = True,
) -> QuantumCircuit:
    """Generate a CNOT gate between data qubit control_qubit and data qubit target_qubit using Bell pairs.
 
    Post processing is used to enable the CNOT gate via the provided classical registers c1 and c2
 
    Assumes that the long-range CNOT gate will be spanning a 1D chain of n-qubits subject to nearest-neighbor
    connections only with the chain starting at the control qubit and finishing at the target qubit.
 
    Assumes that control_qubit < target_qubit (as integers) and that the provided circuit qc has |0> set
    qubits control_qubit+1, ..., target_qubit-1
 
    n = target_qubit - control_qubit - 1 : Number of qubits between the target and control qubits
    k = int(n/2) : Number of Bell pairs created
 
    Args:
        qc (QuantumCicruit) : A Quantum Circuit to add the long range localized unitary CNOT
        control_qubit (int) : The qubit used as the control.
        target_qubi (int) : The qubit targeted by the gate.
 
    Optional Args:
        c1 (ClassicalRegister) : Default = None. Required if n > 1. Register requires k bits
        c2 (ClassicalRegister) : Default = None. Required if n > 0. Register requires n - k bits
        add_barriers (bool) : Default = True. Include barriers before and after long range CNOT
 
    Returns:
        QuantumCircuit
    """
    assert target_qubit > control_qubit
    n = target_qubit - control_qubit - 1
    k = int(n / 2)
 
    # Determine where to start the bell pairs and
    # add an extra CNOT when n is odd
    if n % 2 == 0:
        x0 = 1
    else:
        x0 = 2
        qc.cx(0, 1)
 
    # Create k Bell pairs
    for i in range(k):
        qc.h(x0 + 2 * i)
        qc.cx(x0 + 2 * i, x0 + 2 * i + 1)
 
    # Entangle Bell pairs and data qubits and measure
    for i in range(k + 1):
        qc.cx(x0 - 1 + 2 * i, x0 + 2 * i)
 
    for i in range(1, k + x0):
        qc.h(2 * i + 1 - x0)
        qc.measure(2 * i + 1 - x0, c2[i - 1])
 
    for i in range(k):
        qc.measure(2 * i + x0, c1[i])
 
    if add_barriers is True:
        qc.barrier()
    return qc

Again, utilize the methods prep_P_ij_conj and meas_P_kl to prepare the circuits for Monte Carlo state certification.

def build_circuits_postproc(
    n: int, samples: List[int]
) -> List[QuantumCircuit]:
    """
    Args:
        n (int): Number of qubits between the control and target qubits
    """
    assert n >= 0, "Error: n needs to be a non-negative integer"
    circuits_all = []
 
    qr = QuantumRegister(
        n + 2, name="q"
    )  # Circuit with n qubits between control and target
    cr = ClassicalRegister(
        2, name="cr"
    )  # Classical register for measuring long range CNOT
 
    k = int(n / 2)  # Number of Bell States to be used
    c1 = ClassicalRegister(
        k, name="c1"
    )  # Classical register needed for post processing
    c2 = ClassicalRegister(
        n - k, name="c2"
    )  # Classical register needed for post processing
 
    # 16 Paulis with non-zero expectation value to prepare and measure
    P_lkji = PauliList(
        [
            "IIII",
            "XIXI",
            "IZIZ",
            "XZXZ",
            "YZYI",
            "ZZZI",
            "YIYZ",
            "ZIZZ",
            "XXIX",
            "IXXX",
            "XYIY",
            "IYXY",
            "ZYYX",
            "YYZX",
            "ZXYY",
            "YXZY",
        ]
    )
 
    for sample in samples:
        P_prep = P_lkji[sample][0:2]
        P_meas = P_lkji[sample][2:4]
        if n > 1:
            circuits = [
                QuantumCircuit(qr, cr, c1, c2, name="CNOT") for i in range(4)
            ]
        elif n == 1:
            circuits = [
                QuantumCircuit(qr, cr, c2, name="CNOT") for i in range(4)
            ]
        elif n == 0:
            circuits = [QuantumCircuit(qr, cr, name="CNOT") for i in range(4)]
        circuits = prep_P_ij_conj(
            circuits, P_prep
        )  # Prepare control and target qubits
        # in eigenstates of P_i^* and P_j^* respectively
        circuits = [
            CNOT_postproc(
                qc=circuit, control_qubit=0, target_qubit=n + 1, c1=c1, c2=c2
            )
            for circuit in circuits
        ]  # Add long range CNOT
        circuits = meas_P_kl(
            circuits, P_meas
        )  # Prepare circuits to measure the control and target
        # qubits in P_k and P_l bases respectively
        circuits_all += circuits
    return circuits_all

This is an example of a circuit (including the state certification step):

# Example Circuit
n = 6
sample = [15]
test_post_proc_circuits = build_circuits_postproc(n, sample)
test_post_proc_circuits[1].draw("mpl")

Output:

Output of the previous code cell

Long-range measurement-based CNOT with feedforward

Finally, examine the case where a long-range CNOT gate is implemented using measurement-based CNOT with feedforward (dynamic circuits). In the following figure, on the left is a circuit for a long-range CNOT gate spanning a 1D chain of n-qubits subject to nearest-neighbor connections only. On the right is an equivalent implementable with local CNOT gates, measurement-based CNOT with feedforward (dynamic circuits).

Two equivalent circuits, one using feedforward

The circuit on the right can be implemented as follows:

def CNOT_dyn(
    qc: QuantumCircuit,
    control_qubit: int,
    target_qubit: int,
    c1: Optional[ClassicalRegister] = None,
    c2: Optional[ClassicalRegister] = None,
    add_barriers: Optional[bool] = True,
) -> QuantumCircuit:
    """Generate a CNOT gate between data qubit control_qubit and data qubit target_qubit using Bell Pairs.
 
    Post processing is used to enable the CNOT gate via the provided classical registers c1 and c2
 
    Assumes that the long-range CNOT gate will be spanning a 1D chain of n-qubits subject to nearest-neighbor
    connections only with the chain starting at the control qubit and finishing at the target qubit.
 
    Assumes that control_qubit < target_qubit (as integers) and that the provided circuit qc has |0> set
    qubits control_qubit+1, ..., target_qubit-1
 
    n = target_qubit - control_qubit - 1 : Number of qubits between the target and control qubits
    k = int(n/2) : Number of Bell pairs created
 
    Args:
        qc (QuantumCicruit) : A Quantum Circuit to add the long range localized unitary CNOT
        control_qubit (int) : The qubit used as the control.
        target_qubi (int) : The qubit targeted by the gate.
 
    Optional Args:
        c1 (ClassicalRegister) : Default = None. Required if n > 1. Register requires k bits
        c2 (ClassicalRegister) : Default = None. Required if n > 0. Register requires n - k bits
        add_barriers (bool) : Default = True. Include barriers before and after long range CNOT
 
    Note: This approach uses two if_test statements. A better (more performant) approach is
    to have the parity values combined into a single classical register and then use a switch
    statement. This was done in the associated paper my modifying the qasm file directly. The ability
    to use a switch statement via Qiskit in this way is a future release capability.
 
    Returns:
        QuantumCircuit
    """
    assert target_qubit > control_qubit
    n = target_qubit - control_qubit - 1
    t = int(n / 2)
 
    if add_barriers is True:
        qc.barrier()
 
    # Determine where to start the bell pairs and
    # add an extra CNOT when n is odd
    if n % 2 == 0:
        x0 = 1
    else:
        x0 = 2
        qc.cx(0, 1)
 
    # Create t Bell pairs
    for i in range(t):
        qc.h(x0 + 2 * i)
        qc.cx(x0 + 2 * i, x0 + 2 * i + 1)
 
    # Entangle Bell pairs and data qubits and measure
    for i in range(t + 1):
        qc.cx(x0 - 1 + 2 * i, x0 + 2 * i)
 
    for i in range(1, t + x0):
        if i == 1:
            qc.h(2 * i + 1 - x0)
            qc.measure(2 * i + 1 - x0, c2[i - 1])
            parity_control = expr.lift(c2[i - 1])
        else:
            qc.h(2 * i + 1 - x0)
            qc.measure(2 * i + 1 - x0, c2[i - 1])
            parity_control = expr.bit_xor(c2[i - 1], parity_control)
 
    for i in range(t):
        if i == 0:
            qc.measure(2 * i + x0, c1[i])
            parity_target = expr.lift(c1[i])
        else:
            qc.measure(2 * i + x0, c1[i])
            parity_target = expr.bit_xor(c1[i], parity_target)
 
    if n > 0:
        with qc.if_test(parity_control):
            qc.z(0)
 
    if n > 1:
        with qc.if_test(parity_target):
            qc.x(-1)
 
    if add_barriers is True:
        qc.barrier()
    return qc

Put it together with the Monte Carlo state certification methods prep_P_ij_conj and meas_P_kl:

def build_circuits_dyn(n: int, samples: List[int]) -> List[QuantumCircuit]:
    """ """
    assert n >= 0, "Error: n needs to be a non-negative integer"
    circuits_all = []
 
    qr = QuantumRegister(
        n + 2, name="q"
    )  # Circuit with n qubits between control and target
    cr = ClassicalRegister(
        2, name="cr"
    )  # Classical register for measuring long range CNOT
 
    k = int(n / 2)  # Number of Bell States to be used
    c1 = ClassicalRegister(
        k, name="c1"
    )  # Classical register needed for post processing
    c2 = ClassicalRegister(
        n - k, name="c2"
    )  # Classical register needed for post processing
 
    # 16 Paulis with non-zero expectation value to prepare and measure
    P_lkji = PauliList(
        [
            "IIII",
            "XIXI",
            "IZIZ",
            "XZXZ",
            "YZYI",
            "ZZZI",
            "YIYZ",
            "ZIZZ",
            "XXIX",
            "IXXX",
            "XYIY",
            "IYXY",
            "ZYYX",
            "YYZX",
            "ZXYY",
            "YXZY",
        ]
    )
 
    for sample in samples:
        P_prep = P_lkji[sample][0:2]
        P_meas = P_lkji[sample][2:4]
        if n > 1:
            circuits = [
                QuantumCircuit(qr, cr, c1, c2, name="CNOT") for i in range(4)
            ]
        elif n == 1:
            circuits = [
                QuantumCircuit(qr, cr, c2, name="CNOT") for i in range(4)
            ]
        elif n == 0:
            circuits = [QuantumCircuit(qr, cr, name="CNOT") for i in range(4)]
        circuits = prep_P_ij_conj(
            circuits, P_prep
        )  # Prepare control and target qubits
        # in eigenstates of P_i^* and P_j^* respectively
        circuits = [
            CNOT_dyn(
                qc=circuit, control_qubit=0, target_qubit=n + 1, c1=c1, c2=c2
            )
            for circuit in circuits
        ]  # Add long range CNOT
        circuits = meas_P_kl(
            circuits, P_meas
        )  # Prepare circuits to measure the control and target
        # qubits in P_k and P_l bases respectively
        circuits_all += circuits
    return circuits_all

This results in the following example:

# View an example circuit with Monte Carlo Prep
 
n_qubits = 4
sample = range(16)
example_post_proc_circuits = build_circuits_dyn(n_qubits, sample)
example_post_proc_circuits[16].draw("mpl")

Output:

Output of the previous code cell

Step 2: Optimize problem for quantum hardware execution

Because you have already specified the physical qubit layout and built the circuits with a line topology in mind, there is no need to further optimize the circuits.

Step 3: Execute using Qiskit primitives

In this step you execute the experiment on the specified backend. A lightweight transpilation is done before submission to ensure that the circuits have all the physical parameters needed for execution on the device. This is done in the submit_circuits function, which also splits the circuits into smaller batches to be submitted to the device.

def submit_circuits(
    min_qubits: int,
    max_qubits: int,
    num_circuits_per_job: int,
    qubit_line: List[int],
    coupling_map: Union[CouplingMap, List],
    samples: List[int],
    optimization_level: int,
    backend: Backend,
    shots: int,
    build_circuits: Callable,
    transpile_dynamic: Optional[bool] = True,
    use_dynamic_decoupling: Optional[bool] = True,
    dd_sequence: Optional[List[Gate]] = [XGate(), XGate()],
    durations: Optional[InstructionDurations] = None,
) -> List[str]:
    """
    Submit circuits in appropriate batches
    """
 
    # Calculated constants and storage variables
    line_length = len(qubit_line)
    num_samples = len(samples)
    num_circuits = (max_qubits - min_qubits + 1) * 4 * num_samples
    nr_jobs = int(num_circuits / num_circuits_per_job)
 
    # Run some parameter checks
    # Min number of qubits between control and target must be a non-negative integer
    assert min_qubits >= 0, "Error: min_qubits must be >= 0"
 
    # Max number of qubits between control and target musts be <= line_length - 2
    assert (
        max_qubits + 2
    ) <= line_length, "Error: max_qubits must be <= len(qubit_line) - 2"
 
    # (max_qubits - min_qubits) must equal to 3(mod 4)
    rem = (max_qubits - min_qubits) % 4
    assert rem == 3, "Fail: (max_qubits - min_qubits) must equal to 3(mod 4)"
 
    # First transpile all the circuits
    print("Transpiling circuits...")
 
    all_transpiled_circs = []
 
    for n in range(min_qubits, max_qubits + 1):
        layout = qubit_line[: n + 2]
        circuits = build_circuits(n, samples)
 
        clear_output(wait=True)
        percentage_completed = (n - min_qubits + 1) / (
            max_qubits - min_qubits + 1
        )
 
        print(
            f"[{percentage_completed:.0%} completed] Transpiling circuits "
            + f"with {n} qubits between CNOT"
        )
 
        # Generate the main Qiskit transpile passes.
        pm = generate_preset_pass_manager(
            coupling_map=coupling_map,
            initial_layout=layout,
            optimization_level=optimization_level,
            backend=backend,
        )
 
        if use_dynamic_decoupling is True:
            # Configure the as-late-as-possible scheduling pass and DD insertion pass
            pm.scheduling = PassManager(
                [
                    ALAPScheduleAnalysis(durations),
                    PadDynamicalDecoupling(durations, dd_sequence),
                ]
            )
 
        transpiled_circuits = pm.run(circuits)
        all_transpiled_circs.extend(transpiled_circuits)
 
    clear_output(wait=True)
    print("Sumbitting jobs ...")
 
    job_ids = []
 
    with Batch(backend=backend):
        sampler = Sampler()
        for job_num in range(nr_jobs):
            transpiled_circs = all_transpiled_circs[
                num_circuits_per_job * job_num : num_circuits_per_job
                * (job_num + 1)
            ]
 
            # Submit circuits
            print("Submitting circuits:")
 
            percentage_completed = job_num / nr_jobs
            print(f"[{percentage_completed:.0%} completed]")
 
            job = sampler.run(transpiled_circs, shots=shots)
            job_ids.append(job.job_id())
            print(
                "Job id for circuits "
                + f"[{num_circuits_per_job*nr_jobs}, {num_circuits_per_job*(nr_jobs + 1) -1 }] : {job.job_id()}"
            )
 
            clear_output(wait=True)
 
        clear_output(wait=True)
    print("All jobs submitted.\n")
 
    # Display qubit ranges and job ids
    for job_num in range(nr_jobs):
        print(
            f"[{num_circuits_per_job*job_num}, {num_circuits_per_job*(job_num + 1)}]: "
            f"Id = {job_ids[job_num]}"
        )
 
    return job_ids

First, set the parameters for the unitary approach and submit circuits.

# Set local parameters
SAMPLES_UNI = SAMPLES
OPTIMIZATION_LEVEL_UNI = OPTIMIZATION_LEVEL
SHOTS_UNI = SHOTS
MIN_NUMBER_QUBITS_UNI = MIN_NUMBER_QUBITS
MAX_NUMBER_QUBITS_UNI = MAX_NUMBER_QUBITS
NUM_CIRCUITS_PER_JOB_UNI = 256
USE_DYNAMIC_DECOUPLING_UNI = False
 
# Submit jobs for using unitary circuit approach
job_ids_uni = submit_circuits(
    MIN_NUMBER_QUBITS_UNI,
    MAX_NUMBER_QUBITS_UNI,
    NUM_CIRCUITS_PER_JOB_UNI,
    QUBIT_LINE,
    COUPLING_MAP_1D,
    SAMPLES_UNI,
    OPTIMIZATION_LEVEL_UNI,
    backend,
    SHOTS_UNI,
    build_circuits_uni,
    use_dynamic_decoupling=USE_DYNAMIC_DECOUPLING_UNI,
)

Output:

All jobs submitted.

[0, 256]: Id = d03qhcxnhqag008vm5m0
[256, 512]: Id = d03qhdxrxz8g008g1bng
[512, 768]: Id = d03qhexqnmvg0082cnz0
[768, 1024]: Id = d03qhg6rxz8g008g1bp0
[1024, 1280]: Id = d03qhhed8drg008jr4bg
[1280, 1536]: Id = d03qhk6qnmvg0082cnzg
[1536, 1792]: Id = d03qhmyd8drg008jr4dg
[1792, 2048]: Id = d03qhqynhqag008vm5n0
[2048, 2304]: Id = d03qhszkzhn0008e4njg
[2304, 2560]: Id = d03qhvqkzhn0008e4nk0
[2560, 2816]: Id = d03qhyf6rr3g008s8s00
[2816, 3072]: Id = d03qj10d8drg008jr4fg
[3072, 3328]: Id = d03qj3rrxz8g008g1br0
[3328, 3584]: Id = d03qj6g6rr3g008s8s1g
[3584, 3840]: Id = d03qj9hd8drg008jr4h0
[3840, 4096]: Id = d03qjcsnhqag008vm5pg

Then, do the same for the measurement-based post-selection approach.

# Set local parameters
SAMPLES_POSTPROC = SAMPLES
OPTIMIZATION_LEVEL_POSTPROC = OPTIMIZATION_LEVEL
SHOTS_POSTPROC = SHOTS
MIN_NUMBER_QUBITS_POSTPROC = MIN_NUMBER_QUBITS
MAX_NUMBER_QUBITS_POSTPROC = MAX_NUMBER_QUBITS
NUM_CIRCUITS_PER_JOB_POSTPROC = 128
USE_DYNAMIC_DECOUPLING_POSTPROC = USE_DYNAMIC_DECOUPLING
DURATIONS_POSTPROC = DURATIONS
 
# Submit jobs for the measurement based post selection approach
job_ids_postproc = submit_circuits(
    MIN_NUMBER_QUBITS_POSTPROC,
    MAX_NUMBER_QUBITS_POSTPROC,
    NUM_CIRCUITS_PER_JOB_POSTPROC,
    QUBIT_LINE,
    COUPLING_MAP_1D,
    SAMPLES_POSTPROC,
    OPTIMIZATION_LEVEL_POSTPROC,
    backend,
    SHOTS_POSTPROC,
    build_circuits_postproc,
    use_dynamic_decoupling=USE_DYNAMIC_DECOUPLING_POSTPROC,
    durations=DURATIONS_POSTPROC,
)

Output:

All jobs submitted.

[0, 128]: Id = d03qm4gnhqag008vm5wg
[128, 256]: Id = d03qm58nhqag008vm5xg
[256, 384]: Id = d03qm5r6rr3g008s8sa0
[384, 512]: Id = d03qm6gqnmvg0082cp8g
[512, 640]: Id = d03qm78d8drg008jr4p0
[640, 768]: Id = d03qm81qnmvg0082cp90
[768, 896]: Id = d03qm8hd8drg008jr4q0
[896, 1024]: Id = d03qm99kzhn0008e4ns0
[1024, 1152]: Id = d03qma1qnmvg0082cp9g
[1152, 1280]: Id = d03qmasrxz8g008g1c0g
[1280, 1408]: Id = d03qmbh6rr3g008s8sag
[1408, 1536]: Id = d03qmd9rxz8g008g1c10
[1536, 1664]: Id = d03qme9rxz8g008g1c20
[1664, 1792]: Id = d03qmf1kzhn0008e4ntg
[1792, 1920]: Id = d03qmgaqnmvg0082cpag
[1920, 2048]: Id = d03qmh2qnmvg0082cpb0
[2048, 2176]: Id = d03qmj2qnmvg0082cpbg
[2176, 2304]: Id = d03qmjtqnmvg0082cpc0
[2304, 2432]: Id = d03qmkjkzhn0008e4nvg
[2432, 2560]: Id = d03qmmt6rr3g008s8sc0
[2560, 2688]: Id = d03qmntnhqag008vm5zg
[2688, 2816]: Id = d03qmptnhqag008vm600
[2816, 2944]: Id = d03qmqtkzhn0008e4nwg
[2944, 3072]: Id = d03qmrvnhqag008vm610
[3072, 3200]: Id = d03qmsvkzhn0008e4nx0
[3200, 3328]: Id = d03qmtvd8drg008jr4sg
[3328, 3456]: Id = d03qmyvkzhn0008e4ny0
[3456, 3584]: Id = d03qmzvkzhn0008e4nyg
[3584, 3712]: Id = d03qn14rxz8g008g1c80
[3712, 3840]: Id = d03qn24qnmvg0082cpd0
[3840, 3968]: Id = d03qn3md8drg008jr4tg
[3968, 4096]: Id = d03qn4wkzhn0008e4nzg

Finally, for the measurement-based dynamic circuit approach:

# Set local parameters
SAMPLES_DYN = SAMPLES
OPTIMIZATION_LEVEL_DYN = OPTIMIZATION_LEVEL
SHOTS_DYN = SHOTS
MIN_NUMBER_QUBITS_DYN = MIN_NUMBER_QUBITS
MAX_NUMBER_QUBITS_DYN = MAX_NUMBER_QUBITS
DURATIONS_DYN = DURATIONS
DD_SEQUENCE_DYN = DD_SEQUENCE
NUM_CIRCUITS_PER_JOB_DYN = 16
USE_DYNAMIC_DECOUPLING_DYN = USE_DYNAMIC_DECOUPLING
 
# Submit jobs for the measurement based dynamic circuit approach
job_ids_dyn = submit_circuits(
    MIN_NUMBER_QUBITS_DYN,
    MAX_NUMBER_QUBITS_DYN,
    NUM_CIRCUITS_PER_JOB_DYN,
    QUBIT_LINE,
    COUPLING_MAP_1D,
    SAMPLES_DYN,
    OPTIMIZATION_LEVEL_DYN,
    backend,
    SHOTS_DYN,
    build_circuits_dyn,
    durations=DURATIONS_DYN,
)

Output:

All jobs submitted.

[0, 16]: Id = d03qqrqkzhn0008e4p7g
[16, 32]: Id = d03qqs7d8drg008jr540
[32, 48]: Id = d03qqsqqnmvg0082cpng
[48, 64]: Id = d03qqt7qnmvg0082cpp0
[64, 80]: Id = d03qqtqd8drg008jr54g
[80, 96]: Id = d03qqtzd8drg008jr550
[96, 112]: Id = d03qqvfrxz8g008g1chg
[112, 128]: Id = d03qqvzqnmvg0082cppg
[128, 144]: Id = d03qqwf6rr3g008s8sp0
[144, 160]: Id = d03qqwznhqag008vm6a0
[160, 176]: Id = d03qqx7kzhn0008e4p8g
[176, 192]: Id = d03qqxqqnmvg0082cpq0
[192, 208]: Id = d03qqy7d8drg008jr55g
[208, 224]: Id = d03qqyqnhqag008vm6ag
[224, 240]: Id = d03qqyznhqag008vm6b0
[240, 256]: Id = d03qqzfnhqag008vm6bg
[256, 272]: Id = d03qqzqkzhn0008e4p90
[272, 288]: Id = d03qr00d8drg008jr560
[288, 304]: Id = d03qr0grxz8g008g1cj0
[304, 320]: Id = d03qr10nhqag008vm6c0
[320, 336]: Id = d03qr1grxz8g008g1cjg
[336, 352]: Id = d03qr1rd8drg008jr570
[352, 368]: Id = d03qr28qnmvg0082cprg
[368, 384]: Id = d03qr2rrxz8g008g1ckg
[384, 400]: Id = d03qr30kzhn0008e4pa0
[400, 416]: Id = d03qr3g6rr3g008s8sq0
[416, 432]: Id = d03qr40qnmvg0082cps0
[432, 448]: Id = d03qr48d8drg008jr57g
[448, 464]: Id = d03qr4rqnmvg0082cpsg
[464, 480]: Id = d03qr58d8drg008jr580
[480, 496]: Id = d03qr5grxz8g008g1cm0
[496, 512]: Id = d03qr60kzhn0008e4pag
[512, 528]: Id = d03qr6g6rr3g008s8sr0
[528, 544]: Id = d03qr70d8drg008jr58g
[544, 560]: Id = d03qr786rr3g008s8srg
[560, 576]: Id = d03qr7rnhqag008vm6dg
[576, 592]: Id = d03qr89nhqag008vm6e0
[592, 608]: Id = d03qr8sqnmvg0082cptg
[608, 624]: Id = d03qr91qnmvg0082cpv0
[624, 640]: Id = d03qr9hnhqag008vm6eg
[640, 656]: Id = d03qra1nhqag008vm6f0
[656, 672]: Id = d03qrahqnmvg0082cpvg
[672, 688]: Id = d03qraskzhn0008e4pb0
[688, 704]: Id = d03qrb9qnmvg0082cpw0
[704, 720]: Id = d03qrbsd8drg008jr5a0
[720, 736]: Id = d03qrc96rr3g008s8ss0
[736, 752]: Id = d03qrcsqnmvg0082cpwg
[752, 768]: Id = d03qrd9d8drg008jr5ag
[768, 784]: Id = d03qrdskzhn0008e4pbg
[784, 800]: Id = d03qre1d8drg008jr5b0
[800, 816]: Id = d03qrehqnmvg0082cpx0
[816, 832]: Id = d03qrf1kzhn0008e4pc0
[832, 848]: Id = d03qrf9rxz8g008g1cn0
[848, 864]: Id = d03qrfs6rr3g008s8st0
[864, 880]: Id = d03qrgad8drg008jr5bg
[880, 896]: Id = d03qrgtqnmvg0082cpxg
[896, 912]: Id = d03qrh2d8drg008jr5c0
[912, 928]: Id = d03qrhjqnmvg0082cpy0
[928, 944]: Id = d03qrj2d8drg008jr5cg
[944, 960]: Id = d03qrjjd8drg008jr5d0
[960, 976]: Id = d03qrjtkzhn0008e4pcg
[976, 992]: Id = d03qrkakzhn0008e4pd0
[992, 1008]: Id = d03qrktd8drg008jr5dg
[1008, 1024]: Id = d03qrm2d8drg008jr5e0
[1024, 1040]: Id = d03qrmjnhqag008vm6h0
[1040, 1056]: Id = d03qrn2qnmvg0082cpzg
[1056, 1072]: Id = d03qrnj6rr3g008s8svg
[1072, 1088]: Id = d03qrntkzhn0008e4pdg
[1088, 1104]: Id = d03qrpanhqag008vm6hg
[1104, 1120]: Id = d03qrpt6rr3g008s8swg
[1120, 1136]: Id = d03qrqanhqag008vm6j0
[1136, 1152]: Id = d03qrqtqnmvg0082cq00
[1152, 1168]: Id = d03qrrbrxz8g008g1cpg
[1168, 1184]: Id = d03qrrvrxz8g008g1cq0
[1184, 1200]: Id = d03qrs3kzhn0008e4pe0
[1200, 1216]: Id = d03qrskqnmvg0082cq0g
[1216, 1232]: Id = d03qrt3qnmvg0082cq10
[1232, 1248]: Id = d03qrtknhqag008vm6jg
[1248, 1264]: Id = d03qrv3rxz8g008g1cqg
[1264, 1280]: Id = d03qrvb6rr3g008s8sx0
[1280, 1296]: Id = d03qrvv6rr3g008s8sxg
[1296, 1312]: Id = d03qrwbkzhn0008e4peg
[1312, 1328]: Id = d03qrwvqnmvg0082cq20
[1328, 1344]: Id = d03qrxkrxz8g008g1cs0
[1344, 1360]: Id = d03qrxvqnmvg0082cq2g
[1360, 1376]: Id = d03qrybqnmvg0082cq30
[1376, 1392]: Id = d03qryvkzhn0008e4pgg
[1392, 1408]: Id = d03qrzbd8drg008jr5fg
[1408, 1424]: Id = d03qrzvqnmvg0082cq3g
[1424, 1440]: Id = d03qs0ckzhn0008e4ph0
[1440, 1456]: Id = d03qs0w6rr3g008s8sz0
[1456, 1472]: Id = d03qs1cnhqag008vm6m0
[1472, 1488]: Id = d03qs1mkzhn0008e4phg
[1488, 1504]: Id = d03qs24kzhn0008e4pj0
[1504, 1520]: Id = d03qs2mnhqag008vm6mg
[1520, 1536]: Id = d03qs2wrxz8g008g1csg
[1536, 1552]: Id = d03qs3crxz8g008g1ct0
[1552, 1568]: Id = d03qs44nhqag008vm6n0
[1568, 1584]: Id = d03qs4md8drg008jr5g0
[1584, 1600]: Id = d03qs4wqnmvg0082cq40
[1600, 1616]: Id = d03qs5mqnmvg0082cq4g
[1616, 1632]: Id = d03qs64rxz8g008g1cv0
[1632, 1648]: Id = d03qs6cd8drg008jr5gg
[1648, 1664]: Id = d03qs74nhqag008vm6p0
[1664, 1680]: Id = d03qs7md8drg008jr5hg
[1680, 1696]: Id = d03qs85d8drg008jr5j0
[1696, 1712]: Id = d03qs8dd8drg008jr5jg
[1712, 1728]: Id = d03qs8xqnmvg0082cq5g
[1728, 1744]: Id = d03qs9drxz8g008g1cvg
[1744, 1760]: Id = d03qs9xrxz8g008g1cw0
[1760, 1776]: Id = d03qsa5nhqag008vm6qg
[1776, 1792]: Id = d03qsan6rr3g008s8t00
[1792, 1808]: Id = d03qsb56rr3g008s8t0g
[1808, 1824]: Id = d03qsbn6rr3g008s8t10
[1824, 1840]: Id = d03qsc5rxz8g008g1cwg
[1840, 1856]: Id = d03qscnd8drg008jr5k0
[1856, 1872]: Id = d03qsd56rr3g008s8t1g
[1872, 1888]: Id = d03qsdnd8drg008jr5kg
[1888, 1904]: Id = d03qse5nhqag008vm6rg
[1904, 1920]: Id = d03qsed6rr3g008s8t20
[1920, 1936]: Id = d03qsexd8drg008jr5m0
[1936, 1952]: Id = d03qsfdd8drg008jr5mg
[1952, 1968]: Id = d03qsfxnhqag008vm6s0
[1968, 1984]: Id = d03qsg6kzhn0008e4pm0
[1984, 2000]: Id = d03qsgp6rr3g008s8t2g
[2000, 2016]: Id = d03qsh6qnmvg0082cq80
[2016, 2032]: Id = d03qshpqnmvg0082cq8g
[2032, 2048]: Id = d03qsj6qnmvg0082cq90
[2048, 2064]: Id = d03qsjpd8drg008jr5ng
[2064, 2080]: Id = d03qsk6rxz8g008g1cxg
[2080, 2096]: Id = d03qskekzhn0008e4pmg
[2096, 2112]: Id = d03qsky6rr3g008s8t3g
[2112, 2128]: Id = d03qsmerxz8g008g1cy0
[2128, 2144]: Id = d03qsmyd8drg008jr5p0
[2144, 2160]: Id = d03qsnenhqag008vm6t0
[2160, 2176]: Id = d03qsnyrxz8g008g1cyg
[2176, 2192]: Id = d03qsped8drg008jr5pg
[2192, 2208]: Id = d03qspyd8drg008jr5q0
[2208, 2224]: Id = d03qsqerxz8g008g1cz0
[2224, 2240]: Id = d03qsqy6rr3g008s8t4g
[2240, 2256]: Id = d03qsrf6rr3g008s8t50
[2256, 2272]: Id = d03qsrzqnmvg0082cq9g
[2272, 2288]: Id = d03qss7nhqag008vm6tg
[2288, 2304]: Id = d03qssq6rr3g008s8t5g
[2304, 2320]: Id = d03qst76rr3g008s8t60
[2320, 2336]: Id = d03qstqd8drg008jr5qg
[2336, 2352]: Id = d03qsv7d8drg008jr5r0
[2352, 2368]: Id = d03qsvqrxz8g008g1d00
[2368, 2384]: Id = d03qsw7d8drg008jr5rg
[2384, 2400]: Id = d03qswfnhqag008vm6v0
[2400, 2416]: Id = d03qswz6rr3g008s8t7g
[2416, 2432]: Id = d03qsxf6rr3g008s8t80
[2432, 2448]: Id = d03qsxzkzhn0008e4png
[2448, 2464]: Id = d03qsyfqnmvg0082cqb0
[2464, 2480]: Id = d03qsyzqnmvg0082cqbg
[2480, 2496]: Id = d03qszfd8drg008jr5s0
[2496, 2512]: Id = d03qszzrxz8g008g1d0g
[2512, 2528]: Id = d03qt086rr3g008s8t9g
[2528, 2544]: Id = d03qt0r6rr3g008s8ta0
[2544, 2560]: Id = d03qt10qnmvg0082cqc0
[2560, 2576]: Id = d03qt1g6rr3g008s8tb0
[2576, 2592]: Id = d03qt28rxz8g008g1d10
[2592, 2608]: Id = d03qt2rkzhn0008e4pq0
[2608, 2624]: Id = d03qt306rr3g008s8tbg
[2624, 2640]: Id = d03qt3gnhqag008vm6wg
[2640, 2656]: Id = d03qt40d8drg008jr5tg
[2656, 2672]: Id = d03qt4gnhqag008vm6x0
[2672, 2688]: Id = d03qt50nhqag008vm6xg
[2688, 2704]: Id = d03qt5gd8drg008jr5v0
[2704, 2720]: Id = d03qt60rxz8g008g1d2g
[2720, 2736]: Id = d03qt6grxz8g008g1d30
[2736, 2752]: Id = d03qt70kzhn0008e4pr0
[2752, 2768]: Id = d03qt7rrxz8g008g1d3g
[2768, 2784]: Id = d03qt89rxz8g008g1d40
[2784, 2800]: Id = d03qt8skzhn0008e4prg
[2800, 2816]: Id = d03qt91kzhn0008e4ps0
[2816, 2832]: Id = d03qt9hkzhn0008e4pt0
[2832, 2848]: Id = d03qta1kzhn0008e4ptg
[2848, 2864]: Id = d03qta9d8drg008jr5vg
[2864, 2880]: Id = d03qtas6rr3g008s8tc0
[2880, 2896]: Id = d03qtb9nhqag008vm6yg
[2896, 2912]: Id = d03qtbsqnmvg0082cqd0
[2912, 2928]: Id = d03qtch6rr3g008s8tcg
[2928, 2944]: Id = d03qtd1nhqag008vm6z0
[2944, 2960]: Id = d03qtdhqnmvg0082cqdg
[2960, 2976]: Id = d03qte1kzhn0008e4pv0
[2976, 2992]: Id = d03qte9rxz8g008g1d50
[2992, 3008]: Id = d03qtesqnmvg0082cqe0
[3008, 3024]: Id = d03qtf9kzhn0008e4pvg
[3024, 3040]: Id = d03qtfsnhqag008vm6zg
[3040, 3056]: Id = d03qtgaqnmvg0082cqeg
[3056, 3072]: Id = d03qtgt6rr3g008s8td0
[3072, 3088]: Id = d03qthaqnmvg0082cqfg
[3088, 3104]: Id = d03qthtkzhn0008e4pw0
[3104, 3120]: Id = d03qtja6rr3g008s8tdg
[3120, 3136]: Id = d03qtjtqnmvg0082cqg0
[3136, 3152]: Id = d03qtkarxz8g008g1d70
[3152, 3168]: Id = d03qtkt6rr3g008s8teg
[3168, 3184]: Id = d03qtma6rr3g008s8tf0
[3184, 3200]: Id = d03qtmtd8drg008jr5w0
[3200, 3216]: Id = d03qtna6rr3g008s8tfg
[3216, 3232]: Id = d03qtntkzhn0008e4px0
[3232, 3248]: Id = d03qtpanhqag008vm710
[3248, 3264]: Id = d03qtptd8drg008jr5wg
[3264, 3280]: Id = d03qtqad8drg008jr5x0
[3280, 3296]: Id = d03qtqtrxz8g008g1d7g
[3296, 3312]: Id = d03qtrkd8drg008jr5xg
[3312, 3328]: Id = d03qtrvrxz8g008g1d80
[3328, 3344]: Id = d03qtsbkzhn0008e4pxg
[3344, 3360]: Id = d03qtsv6rr3g008s8tg0
[3360, 3376]: Id = d03qttvqnmvg0082cqh0
[3376, 3392]: Id = d03qtvbnhqag008vm71g
[3392, 3408]: Id = d03qtvvrxz8g008g1d8g
[3408, 3424]: Id = d03qtwb6rr3g008s8tgg
[3424, 3440]: Id = d03qtwvkzhn0008e4pz0
[3440, 3456]: Id = d03qtxbd8drg008jr5yg
[3456, 3472]: Id = d03qtxvkzhn0008e4pzg
[3472, 3488]: Id = d03qtybrxz8g008g1d90
[3488, 3504]: Id = d03qtyvnhqag008vm720
[3504, 3520]: Id = d03qtzb6rr3g008s8th0
[3520, 3536]: Id = d03qv04nhqag008vm72g
[3536, 3552]: Id = d03qv0mnhqag008vm730
[3552, 3568]: Id = d03qv14qnmvg0082cqj0
[3568, 3584]: Id = d03qv1mkzhn0008e4q0g
[3584, 3600]: Id = d03qv24nhqag008vm740
[3600, 3616]: Id = d03qv2mnhqag008vm74g
[3616, 3632]: Id = d03qv34nhqag008vm750
[3632, 3648]: Id = d03qv3mkzhn0008e4q1g
[3648, 3664]: Id = d03qv44nhqag008vm75g
[3664, 3680]: Id = d03qv4md8drg008jr5zg
[3680, 3696]: Id = d03qv54rxz8g008g1db0
[3696, 3712]: Id = d03qv5wqnmvg0082cqjg
[3712, 3728]: Id = d03qv6c6rr3g008s8tkg
[3728, 3744]: Id = d03qv6wkzhn0008e4q3g
[3744, 3760]: Id = d03qv7crxz8g008g1dc0
[3760, 3776]: Id = d03qv7wqnmvg0082cqk0
[3776, 3792]: Id = d03qv8nqnmvg0082cqm0
[3792, 3808]: Id = d03qv8xrxz8g008g1dcg
[3808, 3824]: Id = d03qv9nrxz8g008g1dd0
[3824, 3840]: Id = d03qva56rr3g008s8tm0
[3840, 3856]: Id = d03qvanrxz8g008g1ddg
[3856, 3872]: Id = d03qvb5qnmvg0082cqmg
[3872, 3888]: Id = d03qvbxqnmvg0082cqn0
[3888, 3904]: Id = d03qvcdrxz8g008g1de0
[3904, 3920]: Id = d03qvcxd8drg008jr60g
[3920, 3936]: Id = d03qvdnqnmvg0082cqng
[3936, 3952]: Id = d03qve5nhqag008vm77g
[3952, 3968]: Id = d03qvend8drg008jr610
[3968, 3984]: Id = d03qvfdrxz8g008g1deg
[3984, 4000]: Id = d03qvfxkzhn0008e4q50
[4000, 4016]: Id = d03qvge6rr3g008s8tn0
[4016, 4032]: Id = d03qvgyqnmvg0082cqp0
[4032, 4048]: Id = d03qvhe6rr3g008s8tng
[4048, 4064]: Id = d03qvhy6rr3g008s8tp0
[4064, 4080]: Id = d03qvjyrxz8g008g1dfg
[4080, 4096]: Id = d03qvkeqnmvg0082cqq0

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

After the experiments have successfully executed, proceed to post-process the resulting counts to gain insight on the final results. You can take advantage of resampling techniques (also known as bootstrapping) to calculate average fidelities and deviations from the experimental counts.

def resample_single_dictionary(d):
    """Resample a single dictionary based on its weights."""
    keys = list(d.keys())
    weights = list(d.values())
    sum_weights = sum(weights)
 
    if sum_weights == 0:
        return d
 
    resampled_keys = random.choices(keys, weights=weights, k=sum_weights)
 
    # Count the occurrences of each key in the resampled keys
    resampled_counts = {}
    for key in resampled_keys:
        resampled_counts[key] = resampled_counts.get(key, 0) + 1
 
    return resampled_counts
 
 
def resample_dict_list(dict_list, n_samples):
    """Resample the entire list of dictionaries n_samples times."""
    resampled_lists = []
 
    for _ in range(n_samples):
        new_version = [resample_single_dictionary(d) for d in dict_list]
        resampled_lists.append(new_version)
 
    return resampled_lists

In addition, to post-process the results, you need to extract the information from the Monte Carlo state certification protocol - thus, depending on the preparation/measurement basis, you will group the results differently. The utility functions below are meant to carry out this procedure:

def parity(string: str) -> int:
    return string.count("1") % 2
 
 
def parities(string: str) -> str:
    strings = string.split()
    parities = [parity(val) for val in strings]
    return parities
 
 
def postproc_counts(counts, i, samples):
    P_lkji = PauliList(
        [
            "IIII",
            "XIXI",
            "IZIZ",
            "XZXZ",
            "YZYI",
            "ZZZI",
            "YIYZ",
            "ZIZZ",
            "XXIX",
            "IXXX",
            "XYIY",
            "IYXY",
            "ZYYX",
            "YYZX",
            "ZXYY",
            "YXZY",
        ]
    )
 
    PauliI = Pauli("I")
    PauliX = Pauli("X")
    PauliZ = Pauli("Z")
 
    P_k = P_lkji[samples[i]][2]
    P_l = P_lkji[samples[i]][3]
 
    # determine parities
    counts_post = {"00": 0, "01": 0, "10": 0, "11": 0}
 
    for key in counts:
        parities_list = parities(key)
        w = len(parities_list)
        if w == 3:
            parity_of_c2, parity_of_c1, _ = parities_list
        elif w == 2:
            parity_of_c1 = 0
            parity_of_c2, _ = parities_list
        else:
            parity_of_c1 = 0
            parity_of_c2 = 0
 
        # add parity_of_c2 to q0 (key[-1]) only if P_k is 'X' or 'Y'
        if P_k == PauliI or P_k == PauliZ:
            parity_of_c2 = 0
 
        # add parity_c1 to q1 (key[-2]) only if P_l is 'I' or 'Z' or 'Y'
        if P_l == PauliX:
            parity_of_c1 = 0
 
        control_qubit_value = int(key[-1])  # Control qubit q0
        target_qubit_value = int(key[-2])  # Target qubit q1
 
        new_control_qubit_value = (control_qubit_value + parity_of_c2) % 2
        new_target_qubit_value = (target_qubit_value + parity_of_c1) % 2
 
        new_key = str(new_target_qubit_value) + str(new_control_qubit_value)
 
        counts_post[new_key] += counts[key]
 
    return counts_post
 
 
def post_process_postproc(count, i, p, q, samples):
    return postproc_counts(count, i, samples)
 
 
def post_process_dyn(count, i, p, q, samples):
    return marginal_counts(count, indices=range(2))
 
 
def process_fidelities(
    counts: Union[dict[str, int], List[dict[str, int]]],
    samples: List[int],
    shots: int,
    post_process: Optional[Callable] = None,
) -> List[float]:
    """Calculate the estimated process fidelities from experiment counts data
 
    Args:
        counts (dict[str:int] or List[dict[str:int]]): counts data from an experiment
        samples (List[int]): which of the 16 Paulis with non-zero expectation value to prepare and measure
        shots (int): Number of shots used in experiment
        post_process (Callable): Post process the counts with post_proc if given. Default = None
    """
    exp_all = []
    # 16 Paulis with non-zero expectation value to prepare and measure
    P_lkji = PauliList(
        [
            "IIII",
            "XIXI",
            "IZIZ",
            "XZXZ",
            "YZYI",
            "ZZZI",
            "YIYZ",
            "ZIZZ",
            "XXIX",
            "IXXX",
            "XYIY",
            "IYXY",
            "ZYYX",
            "YYZX",
            "ZXYY",
            "YXZY",
        ]
    )
    sign_rho_lkji = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1]
 
    PauliI = Pauli("I")
 
    for i in range(len(samples)):
        P_i = P_lkji[samples[i]][0]
        P_j = P_lkji[samples[i]][1]
        P_k = P_lkji[samples[i]][2]
        P_l = P_lkji[samples[i]][3]
 
        exp = 0
        # initial state p with eig value p_eig prepared
        for p in range(2):
            if P_i == PauliI:
                p_eig = 1
            else:
                p_eig = (-1) ** p
 
            # initial state q with eig value q_eig prepared
            for q in range(2):
                if P_j == PauliI:
                    q_eig = 1
                else:
                    q_eig = (-1) ** q
 
                # post process count if provided
                if post_process is not None:
                    if len(counts) > 0:
                        counts_post = post_process(
                            counts[i * 4 + 2 * p + q], i, p, q, samples
                        )
                else:
                    if len(counts) > 0:
                        counts_post = counts[i * 4 + 2 * p + q]
 
                # measurement projecting to states r with eig value r_eig
                for r in range(2):
                    if P_k == PauliI:
                        r_eig = 1
                    else:
                        r_eig = (-1) ** r
                    for s in range(2):
                        if P_l == PauliI:
                            s_eig = 1
                        else:
                            s_eig = (-1) ** s
 
                        str_r = str(r)
                        str_s = str(s)
 
                        try:
                            exp += (
                                p_eig
                                * q_eig
                                * s_eig
                                * r_eig
                                * counts_post[str_s + str_r]
                                / shots
                                / 4
                                / sign_rho_lkji[samples[i]]
                            )
                        except KeyError:
                            pass
 
        exp_all.append(exp)
    return exp_all
 
 
def get_counts_from_bitarray(instance):
    """
    Extract counts from result data
    """
    for field, value in instance.__dict__.items():
        if isinstance(value, BitArray):
            return value.get_counts()
    return None
 
 
def cal_average_fidelities(
    job_ids: List[str],
    min_qubits: int,
    max_qubits: int,
    samples: List[int],
    shots: int,
    num_circuits_per_job: int,
    post_process: Optional[Callable] = None,
    all_counts: Optional[List[Dict]] = None,
    display: Optional[bool] = True,
    debug: Optional[bool] = False,
    n_bootstrap_sample: Optional[int] = 4,
) -> (List[float], List[float]):
    """
    Calculate the average gate fidelities
    """
 
    proc_fidelities = []
    proc_std = []
    nr_jobs = len(job_ids)
    num_samples = len(samples)
    empty_counts = {"00": 0, "01": 0, "10": 0, "11": 0}
    if all_counts is None:
        counts_flag = False
        all_counts = []
    else:
        counts_flag = True
 
    if debug is True:
        print(f"{nr_jobs} to process")
    if len(all_counts) == 0:
        for j in range(nr_jobs):
            job = service.job(job_ids[j])
 
            if job.status() == "RUNNING":
                raise ValueError(
                    f"Job id : {job_ids[j]} returned status 'RUNNING', wait for all jobs to complete then try again."
                )
            if job.status() != "DONE":
                print(
                    f"Job id : {job_ids[j]} returned status '{job.status()}'; adding empty dictionaries."
                )
                all_counts += [empty_counts] * num_circuits_per_job
                continue
            if display is True:
                print(
                    f"Retrieving job data: {job_ids[j]}: {j} of {nr_jobs-1}"
                )
            result = job.result()
            for i in range(len(result)):
                counts = get_counts_from_bitarray(result[i].data)
                all_counts.append(counts)
 
            if debug is False:
                clear_output(wait=True)
    else:
        print("Using provided all_counts data instead of loading from server")
        print(all_counts)
 
    for n in range(min_qubits, max_qubits + 1):
        if display is True:
            print(
                f"Resampling counts for n = {n}: {max_qubits + 1 - n} remaining"
            )
        counts = all_counts[
            (n - min_qubits) * 4 * num_samples : (n - min_qubits + 1)
            * 4
            * num_samples
        ]
        proc_fid_temp = []
 
        for _ in range(n_bootstrap_sample):
            resample_counts = resample_dict_list(counts, 1)[0]
            sample_fidelities = process_fidelities(
                resample_counts, samples, shots, post_process
            )
            proc_fid_temp.append(np.mean(sample_fidelities))
 
        mean, std = (
            np.mean(np.array(proc_fid_temp)),
            np.std(np.array(proc_fid_temp)),
        )
        proc_fidelities.append(mean)
        proc_std.append(std)
        if debug is False:
            clear_output(wait=True)
 
    if display is True:
        print("Process fidelities:")
        print(["{0:0.3f}".format(i) for i in proc_fidelities])
        print("Process fidelities std:")
        print(["{0:0.3f}".format(i) for i in proc_std])
 
    # Calculate average gate fidelity from the process fidelity
 
    avg_gate_fidelities = []
 
    for i in range(len(proc_fidelities)):
        # Use result of Horodecki et al. to calculate the average gate fidelity
        avg_gate_fidelity = (proc_fidelities[i] * 4 + 1) / 5
        avg_gate_fidelities.append(avg_gate_fidelity)
 
    if display is True:
        print("Average Gate Fidelites")
        print(["{0:0.3f}".format(i) for i in avg_gate_fidelities])
 
    # Calculate average gate fidelity std from the process fidelity std
 
    avg_gate_stds = []
 
    for i in range(len(proc_std)):
        # We scale the std as in the average gate fidelity
        avg_gate_std = (proc_std[i] * 4) / 5
        avg_gate_stds.append(avg_gate_std)
 
    if display is True:
        print("Average Gate Std")
        print(["{0:0.3f}".format(i) for i in avg_gate_stds])
 
    if counts_flag is True:
        return (avg_gate_fidelities, avg_gate_stds, all_counts)
    else:
        return (avg_gate_fidelities, avg_gate_stds)
# No post processing of the counts is required in the Unitary circuits. The average gate fidelities can now be calculated:
avg_gate_fidelities_uni, avg_gate_stds_uni = cal_average_fidelities(
    job_ids_uni,
    MIN_NUMBER_QUBITS_UNI,
    MAX_NUMBER_QUBITS_UNI,
    SAMPLES_UNI,
    SHOTS_UNI,
    NUM_CIRCUITS_PER_JOB_UNI,
)
 
 
avg_gate_fidelities_postproc, avg_gate_stds_postproc = cal_average_fidelities(
    job_ids_postproc,
    MIN_NUMBER_QUBITS_POSTPROC,
    MAX_NUMBER_QUBITS_POSTPROC,
    SAMPLES_POSTPROC,
    SHOTS_POSTPROC,
    NUM_CIRCUITS_PER_JOB_POSTPROC,
    post_process=post_process_postproc,
)
 
 
avg_gate_fidelities_dyn, avg_gate_stds_dyn = cal_average_fidelities(
    job_ids_dyn,
    MIN_NUMBER_QUBITS_DYN,
    MAX_NUMBER_QUBITS_DYN,
    SAMPLES_DYN,
    SHOTS_DYN,
    NUM_CIRCUITS_PER_JOB_DYN,
    post_process=post_process_dyn,
)

Output:

Process fidelities:
['0.867', '0.614', '0.542', '0.120', '0.590', '0.467', '0.495', '0.505', '0.529', '0.499', '0.499', '0.464', '0.444', '0.402', '0.424', '0.330', '0.371', '0.395', '0.383', '0.385', '0.372', '0.318', '0.333', '0.339', '0.347', '0.335', '0.352', '0.216', '0.310', '0.300', '0.233', '0.283', '0.280', '0.208', '0.256', '0.270', '0.171', '0.202', '0.123', '0.209', '0.199', '0.186', '0.148', '0.164', '0.156', '0.156', '0.174', '0.123', '0.190', '0.176', '0.134', '0.147', '0.074', '0.156', '0.159', '0.170', '0.142', '0.150', '0.141', '0.124', '0.157', '0.151', '0.069', '0.116']
Process fidelities std:
['0.002', '0.003', '0.002', '0.003', '0.003', '0.003', '0.004', '0.003', '0.001', '0.002', '0.004', '0.002', '0.002', '0.002', '0.004', '0.005', '0.004', '0.002', '0.003', '0.002', '0.002', '0.005', '0.003', '0.004', '0.004', '0.005', '0.003', '0.003', '0.003', '0.002', '0.003', '0.004', '0.002', '0.003', '0.007', '0.003', '0.005', '0.001', '0.003', '0.000', '0.001', '0.002', '0.002', '0.003', '0.005', '0.002', '0.003', '0.003', '0.002', '0.003', '0.003', '0.002', '0.002', '0.002', '0.001', '0.002', '0.003', '0.002', '0.002', '0.004', '0.004', '0.005', '0.003', '0.002']
Average Gate Fidelites
['0.894', '0.691', '0.633', '0.296', '0.672', '0.574', '0.596', '0.604', '0.624', '0.600', '0.600', '0.571', '0.555', '0.521', '0.539', '0.464', '0.497', '0.516', '0.507', '0.508', '0.497', '0.454', '0.467', '0.471', '0.478', '0.468', '0.481', '0.373', '0.448', '0.440', '0.386', '0.426', '0.424', '0.366', '0.405', '0.416', '0.337', '0.362', '0.299', '0.367', '0.359', '0.349', '0.319', '0.331', '0.325', '0.324', '0.339', '0.298', '0.352', '0.341', '0.307', '0.317', '0.259', '0.325', '0.327', '0.336', '0.314', '0.320', '0.313', '0.299', '0.326', '0.321', '0.255', '0.293']
Average Gate Std
['0.002', '0.002', '0.002', '0.003', '0.002', '0.003', '0.003', '0.003', '0.001', '0.001', '0.003', '0.002', '0.001', '0.002', '0.003', '0.004', '0.003', '0.001', '0.002', '0.002', '0.002', '0.004', '0.003', '0.003', '0.003', '0.004', '0.003', '0.002', '0.002', '0.001', '0.002', '0.004', '0.001', '0.002', '0.005', '0.002', '0.004', '0.001', '0.003', '0.000', '0.001', '0.001', '0.002', '0.002', '0.004', '0.002', '0.002', '0.002', '0.001', '0.002', '0.003', '0.001', '0.002', '0.002', '0.001', '0.002', '0.003', '0.002', '0.002', '0.003', '0.003', '0.004', '0.002', '0.002']

Plot the results

To appreciate the results visually, the cell below plots the estimated gate fidelities measured at varying distance between entangled qubits for the three different methods. In general, the fidelity will decrease with increasing distance. The results show that although the unitary method (using SWAPs to implement a long-range entangling interaction) performs better at short distances, there is a cross-over to a regime where dynamic circuits become a better option. This is true for both the measurement-and-feedforward technique as well as the post-processing one.

fig, ax = plt.subplots()
 
ax.errorbar(
    range(MIN_NUMBER_QUBITS_UNI, MAX_NUMBER_QUBITS_UNI + 1),
    avg_gate_fidelities_uni,
    avg_gate_stds_uni,
    fmt="o-.",
    color="c",
    label="Unitary",
)
ax.errorbar(
    range(MIN_NUMBER_QUBITS_UNI, MAX_NUMBER_QUBITS_UNI + 1),
    avg_gate_fidelities_postproc,
    avg_gate_stds_postproc,
    fmt="o-.",
    color="y",
    label="Post-processing",
)
ax.errorbar(
    range(MIN_NUMBER_QUBITS_UNI, MAX_NUMBER_QUBITS_UNI + 1),
    avg_gate_fidelities_dyn,
    avg_gate_stds_dyn,
    fmt="o-.",
    color="m",
    label="Dynamic",
)
ax.axhline(y=1 / 4, color="g", linestyle="--", label="Random gate")
legend = ax.legend(frameon=True)
for text in legend.get_texts():
    text.set_color("black")  # Set the legend text color to black
legend.get_frame().set_facecolor(
    "white"
)  # Set the legend background color to white
legend.get_frame().set_edgecolor(
    "black"
)  # Optional: set the legend border color to black
ax.set_xlabel("Number of qubits between control and target", color="black")
ax.set_ylabel("Teleported gate fidelity", color="black")
ax.grid(linestyle=":", linewidth=0.6, alpha=0.4, color="gray")
ax.set_ylim((0.2, 1))
ax.set_facecolor("white")  # Set the background color of the axes
fig.patch.set_facecolor("white")  # Set the background color of the figure
 
# Ensure the axis lines and ticks are visible
for spine in ax.spines.values():
    spine.set_visible(True)
    spine.set_color("black")  # Set the color of the axis lines to black
ax.tick_params(axis="x", colors="black")
ax.tick_params(axis="y", colors="black")
 
plt.show()

Output:

Output of the previous code cell

Data from the paper

The results from this tutorial are likely to vary from the results of the paper due to different calibrations and machines used. The code presented above also uses a slightly different method to calculate parities than was used in the paper, as well as some other differences to make the notebook cleaner and more accessible to a wider audience.

Plot from Paper

Plot of teleported gate fidelity against the number of qubits.

Appendix: Calculating the average fidelity

The fidelity [2] of two states ρ\rho and σ\sigma is defined by

F(ρ,σ)=Tr(ρσρ)2\mathcal{F}(\rho,\sigma) = \mathrm{Tr}\left(\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}} \right)^2

If one of ρ\rho or σ\sigma is a pure state then this reduces to F(ρ,σ)=Tr(ρσ)\mathcal{F}(\rho,\sigma)=\mathrm{Tr}(\rho\sigma). Gate fidelity is a tool for comparing how well the implemented quantum channel ξ\xi approximates the desired unitary channel U(ρ)=UρU†\mathcal{U}(\rho) = U\rho{U^\dagger}. Gate fidelity is a function defined on pure states as follows:

Fξ,U(∣ϕ⟩):=F(ξ(∣ϕ⟩⟨ϕ∣),U(∣ϕ⟩⟨ϕ∣))=⟨ϕ∣(U†∘ξ)(∣ϕ⟩⟨ϕ∣)∣ϕ⟩:=FU†∘ξ(∣ϕ⟩).\mathcal{F}_{\xi,\mathcal{U}}(|\phi\rangle) := \mathcal{F}\bigl(\xi(|\phi\rangle\langle\phi|), \mathcal{U}(|\phi\rangle\langle\phi|)\bigl) = \langle\phi|(\mathcal{U}^\dagger\circ\xi)(|\phi\rangle\langle\phi|)|\phi\rangle := \mathcal{F}_{\mathcal{U}^\dagger\circ\xi}(|\phi\rangle).

Here FU†∘ξ\mathcal{F}_{\mathcal{U}^\dagger\circ\xi} can be thought of as measuring how noisy the channel U†∘ξ\mathcal{U}^\dagger\circ\xi is. The average gate fidelity of a channel U†∘ξ\mathcal{U}^\dagger\circ\xi is defined by averaging the gate fidelity via the induced haar measure (the Fubini-Stufy meaure):

Favg(U,ξ):=Favg(U†∘ξ):=∫⟨ϕ∣(U†∘ξ)(∣ϕ⟩⟨ϕ∣)∣ϕ⟩dϕ\mathcal{F}_{avg}(\mathcal{U},\xi):=\mathcal{F}_{avg}(\mathcal{U}^\dagger\circ\xi) := \int\langle\phi|(\mathcal{U}^\dagger\circ\xi)(|\phi\rangle\langle\phi|)|\phi\rangle d\phi

To calculate the average gate fidelity of the channel U†∘ξ\mathcal{U}^\dagger\circ\xi we use a result of Horodecki et al. [3] which relates the average gate fidelity to the entanglement fedilty of a channel. The entanglement fidelity of the channel U†∘ξ\mathcal{U}^\dagger\circ\xi is defined as

Fent(U†∘ξ):=Fent(ρU†∘ξ):=⟨ψ+∣ρU†∘ξ∣ψ+⟩=Tr(U†∘ξ)/d2.\mathcal{F}_{ent}(\mathcal{U}^\dagger\circ\xi) := \mathcal{F}_{ent}(\rho_{\mathcal{U}^\dagger\circ\xi}) :=\langle\psi_+|\rho_{\mathcal{U}^\dagger\circ\xi}|\psi_+\rangle = \mathrm{Tr}(\mathcal{U}^\dagger\circ\xi)/d^2.

where ρU†∘ξ\rho_{\mathcal{U}^\dagger\circ\xi} is the density operator obtained from the channel U†∘ξ\mathcal{U}^\dagger\circ\xi via the Choi-Jamoiłkawski isomorphism

ρU†∘ξ=(I⊗(U†∘ξ))(∣ψ+⟩⟨ψ+∣)\rho_{\mathcal{U}^\dagger\circ\xi} = \bigl(I\otimes(\mathcal{U}^\dagger\circ\xi)\bigr)(|\psi_+\rangle\langle\psi_+|)

and where ∣ψ+⟩|\psi_+\rangle is the maximally entangle state

∣ψ+⟩=1d∑i=0d−1∣i⟩⊗∣i⟩.|\psi_+\rangle = \frac{1}{\sqrt{d}}\sum_{i=0}^{d-1}|i\rangle \otimes |i\rangle.

In our specific situation, where U\mathcal{U} is a unitary channel, the entanglement fidelity of U†∘ξ\mathcal{U}^\dagger\circ\xi can be written in terms of the process fidelity of the two Choi states ρU\rho_\mathcal{U} and ρξ\rho_{\xi} as follows:

Fent(U†∘ξ)=Fproc(ρU,ρξ):=F(ρU,ρξ)\mathcal{F}_{ent}(\mathcal{U}^\dagger\circ\xi) = \mathcal{F}_{proc}(\rho_\mathcal{U}, \rho_{\xi}) := \mathcal{F}(\rho_\mathcal{U}, \rho_{\xi})

and so we see via Proposition 1 of Horodecki et al. [3] that

Favg(U,ξ)=Favg(U†∘ξ)=dFent(U†∘ξ)+1d+1=dF(ρU,ρξ)+1d+1\mathcal{F}_{avg}(\mathcal{U},\xi) = \mathcal{F}_{avg}(\mathcal{U}^\dagger\circ\xi) = \frac{d\mathcal{F}_{ent}(\mathcal{U}^\dagger\circ\xi) + 1}{d+1} = \frac{d\mathcal{F}(\rho_\mathcal{U}, \rho_{\xi}) +1}{d+1}

Calculating the process fidelity between two states can now be achieved via Monte Carlo state certification.

As per [4] a direct implementation of the quantum Monte Carlo state certification would prepare a maximally entangled state ∣ψ+⟩|\psi_+\rangle, apply ξ\xi to half of the system, and then measure random Pauli operators on all qubits. A more practical approach consists of preparing the complex conjugate of random product of eigenstates of local Pauli operators (corresponding to the resulting state after half of the entangled state is measured destructively), applying the transformation ξ\xi to the system, and finally measuring a random Pauli operator on each qubit. This can be seen from the following equality:

Tr[(Pi⊗Pj⊗Pk⊗Pl)(I⊗ξ)(∣ψ+⟩⟨ψ+∣)]=1dTr[(Pk⊗Pl)⋅ξ(Pi∗⊗Pj∗)]\mathrm{Tr}\bigl[(P_i\otimes P_j\otimes P_k\otimes P_l)(I\otimes\xi)(|\psi_+\rangle\langle\psi_+|)\bigr] = \frac{1}{d}\mathrm{Tr}\bigl[(P_k\otimes P_l)\cdot \xi(P_i^*\otimes P_j^*)\bigl]

The following three experiments use the modified and simplified version of Monte Carlo state certification combined with the relations derived above to calculate the average gate fidelity of the channel ξ\xi. For more details see [1] and associated references.

References

[1] Efficient Long-Range Entanglement using Dynamic Circuits, by Elisa Bäumer, Vinay Tripathi, Derek S. Wang, Patrick Rall, Edward H. Chen, Swarnadeep Majumder, Alireza Seif, Zlatko K. Minev. IBM Quantum, (2023). https://arxiv.org/abs/2308.13065

[2] Quantum Computation and Quantum Information, by Nielsen and Chuang, Section 9.2.2, (2010)

[3] General teleportation channel, singlet fraction, and quasidistillation, by M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Rev. A 60, 1888 (1999).

[4] Practical characterization of quantum devices without tomography, by M. P. da Silva, O. Landon-Cardinal, and D. Poulin, Phys. Rev. Lett. 107, 210404 (2011).

Tutorial survey

Please take one minute to provide feedback on this tutorial. Your insights will help us improve our content offerings and user experience.

Link to survey

© IBM Corp. 2024