Distributed Quantum Computing over Classical Channels: Wire Cutting, Gate Cutting, and Quasiprobability Decomposition

Modern monolithic quantum processing units (QPUs) face a severe physical bottleneck. As qubit counts increase on a single die, cryogenic systems encounter thermal budget limits, cross-talk amplifies across tightly packed drive lines, and monolithic wafer yields drop exponentially.

Building warehouse-scale quantum computers will ultimately require distributed architectures. While long-term roadmaps envision coherent optical interconnects and distributed quantum state teleportation via shared Einstein-Podolsky-Rosen (EPR) pairs, quantum transducer fidelities and optical-matter coupling rates currently restrict coherent inter-QPU links.

This raises a central question: Can we execute a large quantum circuit across physically disconnected QPUs using only classical communication channels?

The answer lies in quasiprobability decomposition (QPD), commonly known as circuit knitting (subsuming wire cutting and gate cutting). By trading coherent quantum memory channels for classical communication and classical sampling overhead, we can evaluate non-local quantum circuits across an arbitrary cluster of independent QPUs.

1. The LOCC Barrier and Quasiprobability Decomposition

In distributed quantum systems, operations are categorized under LOCC (Local Operations and Classical Communication). Local unitaries acting purely within partition $A$ ($\mathcal{H}_A$) or partition $B$ ($\mathcal{H}_B$) are resource-free. However, entangling operations between the two partitions cannot be implemented deterministically via LOCC without pre-shared entanglement.

       Monolithic Execution                    Distributed LOCC Execution
   ┌───────────────────────────┐         ┌───────────┐       ┌───────────┐
   │ QPU A ───■───             │         │   QPU A   │       │   QPU B   │
   │          │   (Coherent)   │  ───►   │  ───[M]   │ ~ ~ ~ │   [P]───  │
   │ QPU B ───■───             │         └───────────┘  LOCC └───────────┘
   └───────────────────────────┘                        Bus

According to the Nielsen-Kempe theorem and fundamental entanglement monotonicity, any bipartite unitary channel $\mathcal{U}{AB}(\rho) = U{AB} \rho U_{AB}^\dagger$ that generates entanglement ($E(\mathcal{U}_{AB}(\rho)) > E(\rho)$) strictly cannot be expressed as a completely positive trace-preserving (CPTP) LOCC map:

$$\mathcal{U}_{AB} \notin \text{LOCC}$$

Overcoming Non-Locality with Quasiprobability

To bypass this constraint without physical quantum channels, we expand the non-local superoperator $\mathcal{E}_{AB}$ as an affine combination of separable, locally implementable CPTP operations ${\mathcal{A}_i \otimes \mathcal{B}_i}$:

$$\mathcal{E}{AB}(\rho) = \sum{i=1}^M c_i \, (\mathcal{A}_i \otimes \mathcal{B}_i)(\rho), \quad c_i \in \mathbb{R}$$

Because $\mathcal{E}_{AB}$ is non-local, the coefficients $c_i$ cannot form a valid classical probability distribution; necessarily, some $c_i < 0$. This representation is a quasiprobability distribution.

We define the total sampling 1-norm $\gamma$:

$$\gamma = \sum_{i=1}^M |c_i| \ge 1$$

We can then define a valid classical probability distribution $p_i$ and an associated phase sign $s_i$:

$$p_i = \frac{|c_i|}{\gamma}, \qquad s_i = \text{sgn}(c_i) \in {-1, +1}$$

Rewriting the non-local expectation value of an observable $\hat{O} = \hat{O}_A \otimes \hat{O}_B$:

$$\langle \hat{O} \rangle = \text{Tr}\Big(\hat{O} \, \mathcal{E}{AB}(\rho)\Big) = \gamma \sum{i=1}^M p_i \, s_i \, \text{Tr}\Big((\hat{O}_A \mathcal{A}_i \otimes \hat{O}_B \mathcal{B}_i)(\rho)\Big)$$

$$\langle \hat{O} \rangle = \gamma \, \mathbb{E}{i \sim p_i} \Big[ s_i \cdot \langle \hat{O}_A \rangle{\mathcal{A}i} \cdot \langle \hat{O}_B \rangle{\mathcal{B}_i} \Big]$$

The Sampling Overhead

To estimate $\langle \hat{O} \rangle$ to within additive precision $\epsilon$ with confidence $1 - \delta$, we sample subcircuit configurations $i$ from $p_i$, execute $\mathcal{A}_i$ on QPU $A$ and $\mathcal{B}_i$ on QPU $B$, record measurement outcomes classically, and compute the signed estimator:

$$\hat{X} = \gamma \cdot s_i \cdot m_A \cdot m_B$$

The variance of this unbiased estimator scales quadratically with $\gamma$:

$$\text{Var}[\hat{X}] = \gamma^2 \, \mathbb{E}[m_A^2 m_B^2] - \langle \hat{O} \rangle^2 \le \gamma^2 - \langle \hat{O} \rangle^2$$

By Chebyshev's and Hoeffding's inequalities, the total number of circuit executions (shots) $N_{\text{shots}}$ scales as:

$$N_{\text{shots}} = \mathcal{O}\left( \frac{\gamma^2}{\epsilon^2} \ln\left(\frac{1}{\delta}\right) \right)$$

The quantity $\gamma^2$ is the sampling overhead. For $K$ independent non-local cuts, the total sampling overhead compounds multiplicatively:

$$\gamma_{\text{total}}^2 = \prod_{k=1}^K \gamma_k^2 = \gamma^{2K}$$

2. Gate Cutting vs. Wire Cutting

Distributed quantum circuit execution relies on two primary cutting primitives: Gate Cutting and Wire Cutting.

          GATE CUTTING                              WIRE CUTTING
       (Decomposing U_AB)                     (Decomposing Channel I)

      QPU A ───[A_i]───                      QPU A ───[Measure M_k]
                 │                                        │
           (Classical LOCC)                         (Classical LOCC)
                 │                                        │
      QPU B ───[B_i]───                      QPU B ──────[Prep S_k]───

A. Gate Cutting

Gate cutting targets the non-local entangling operator directly. Consider a parameterized two-qubit rotation $R_{ZZ}(\theta) = \exp\left(-i \frac{\theta}{2} Z \otimes Z\right)$:

$$R_{ZZ}(\theta) = \cos\left(\frac{\theta}{2}\right) (I \otimes I) - i \sin\left(\frac{\theta}{2}\right) (Z \otimes Z)$$

At the superoperator level, decomposing the channel $\mathcal{E}{R{ZZ}}(\rho) = R_{ZZ}(\theta)\rho R_{ZZ}^\dagger(\theta)$ yields a quasiprobability overhead:

$$\gamma(R_{ZZ}(\theta)) = 1 + 2\left|\sin\left(\frac{\theta}{2}\right)\cos\left(\frac{\theta}{2}\right)\right| = 1 + |\sin(\theta)|$$

For a maximally entangling Controlled-$Z$ ($CZ$) or $CNOT$ gate ($\theta = \pi/2$), the optimal single-gate decomposition yields $\gamma = 3$, leading to a shot overhead of $\gamma^2 = 9$.

B. Wire Cutting (State Slicing)

Wire cutting severs a qubit transmission line across time slices between partitions. The identity channel $\mathcal{I}(\rho) = \rho$ transmitting a single-qubit state from partition $A$ to partition $B$ is expanded across an orthonormal operator basis. Using the Pauli basis ${I, X, Y, Z}$:

$$\mathcal{I}(\rho) = \frac{1}{2} \text{Tr}(\rho) I + \frac{1}{2} \text{Tr}(X\rho) X + \frac{1}{2} \text{Tr}(Y\rho) Y + \frac{1}{2} \text{Tr}(Z\rho) Z$$

Each Pauli measurement $\text{Tr}(\sigma \rho)$ on QPU $A$ is mapped to projective eigenstates $|\sigma_\pm\rangle$, and state preparations on QPU $B$ inject the corresponding eigenstates:

$$P_{\sigma, \pm} = \frac{I \pm \sigma}{2}, \qquad \sigma \in {X, Y, Z}$$

Expressing the wire identity map into local measurement-preparation pairs $(M_k, S_k)$:

$$\mathcal{I}(\rho) = \text{Tr}(P_{Z,+} \rho) P_{Z,+} + \text{Tr}(P_{Z,-} \rho) P_{Z,-} + \frac{1}{2} \sum_{\sigma \in {X, Y}} \sum_{s, s' \in {+, -}} s s' \, \text{Tr}(P_{\sigma, s} \rho) P_{\sigma, s'}$$

Summing the absolute values of the expansion coefficients gives the single-wire sampling overhead:

$$\gamma_{\text{wire}} = 1 + 1 + 4 \times \left(\frac{1}{2}\right) + 4 \times \left(\frac{1}{2}\right) = 4$$

A single wire cut incurs a shot multiplier of $\gamma^2 = 16$. However, wire cutting allows arbitrary multi-qubit subcircuits to run before and after the cut boundary, making it structurally flexible for deep circuits with low spatial interconnect density.

3. Distributed Execution Architecture

Executing a cut quantum circuit across distributed QPUs requires coordinated execution and classical aggregation:

                  ┌─────────────────────────────────┐
                  │ Global Quantum Circuit & DAG    │
                  └────────────────┬────────────────┘
                                   │ Min-Cut / ILP Partitioning
                    ┌──────────────┴──────────────┐
                    ▼                             ▼
         ┌─────────────────────┐       ┌─────────────────────┐
         │ QPU A Subcircuits   │       │ QPU B Subcircuits   │
         │ (Measurements M_i)  │       │ (Preparations S_i)  │
         └──────────┬──────────┘       └──────────┬──────────┘
                    │ Execution                   │ Execution
                    ▼                             ▼
         ┌─────────────────────┐       ┌─────────────────────┐
         │ Bitstrings {m_A}    │       │ Bitstrings {m_B}    │
         └──────────┬──────────┘       └──────────┬──────────┘
                    │                             │
                    └──────────────┬──────────────┘
                                   │ Classical LOCC Bus
                                   ▼
                  ┌─────────────────────────────────┐
                  │ Classical Tensor Reconstruction │
                  │   = γ * Σ s_i *  * │
                  └─────────────────────────────────┘

1. DAG Partitioning: The global quantum directed acyclic graph (DAG) is partitioned across available QPUs using hypergraph min-cut algorithms to minimize the cut cost $\sum \log(\gamma_i)$.
2. QPD Variant Expansion: The cuts are expanded into discrete, parameterized subcircuit batches.
3. Parallel QPU Dispatch: Subcircuits are executed asynchronously across independent QPUs over standard network fabrics (e.g., gRPC over TCP/IP).
4. Classical Tensor Contraction: Measurement bitstrings from QPU $A$ and QPU $B$ are collected and recombined using quasiprobability weight vectors to reconstruct unbiased expectation values.

4. Production-Ready Implementation

The following Python program implements a complete simulation of distributed quantum computing via classical communication. It models a monolithic 2-qubit circuit ($R_y \to \text{Wire Cut} \to R_x \to CNOT \to \text{Measurement}$), splits it across two independent QPUs using wire cutting QPD, executes finite-shot Monte Carlo sampling, and reconstructs the exact expectation value.

"""
Distributed Quantum Computing via Classical Communication
Simulation of Wire Cutting and Quasiprobability Decomposition (QPD).
"""

from dataclasses import dataclass
from typing import List, Tuple
import numpy as np


# =====================================================================
# 1. Fundamental Quantum Algebra and Quantum States
# =====================================================================

I2 = np.eye(2, dtype=np.complex128)
PAULI_X = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128)
PAULI_Y = np.array([[0.0, -1.0j], [1.0j, 0.0]], dtype=np.complex128)
PAULI_Z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.complex128)

# Basis Eigenstates
Z_PLUS = np.array([1.0, 0.0], dtype=np.complex128)
Z_MINUS = np.array([0.0, 1.0], dtype=np.complex128)
X_PLUS = (Z_PLUS + Z_MINUS) / np.sqrt(2.0)
X_MINUS = (Z_PLUS - Z_MINUS) / np.sqrt(2.0)
Y_PLUS = (Z_PLUS + 1.0j * Z_MINUS) / np.sqrt(2.0)
Y_MINUS = (Z_PLUS - 1.0j * Z_MINUS) / np.sqrt(2.0)

# Projector Helper
def projector(state: np.ndarray) -> np.ndarray:
    return np.outer(state, np.conjugate(state))


# =====================================================================
# 2. Gate Definitions
# =====================================================================

def rx_gate(theta: float) -> np.ndarray:
    return np.cos(theta / 2.0) * I2 - 1.0j * np.sin(theta / 2.0) * PAULI_X

def ry_gate(theta: float) -> np.ndarray:
    return np.cos(theta / 2.0) * I2 - np.sin(theta / 2.0) * (1.0j * PAULI_Y)

def cnot_gate() -> np.ndarray:
    return np.array(
        [
            [1.0, 0.0, 0.0, 0.0],
            [0.0, 1.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, 1.0],
            [0.0, 0.0, 1.0, 0.0],
        ],
        dtype=np.complex128,
    )


# =====================================================================
# 3. Wire Cutting Quasiprobability Decomposition Frame
# =====================================================================

@dataclass(frozen=True)
class QPDComponent:
    """Represents a single quasi-probability channel element."""
    weight: float
    meas_projector: np.ndarray  # Operator applied to QPU A cut boundary
    prep_state: np.ndarray      # State prepared on QPU B cut boundary


def get_wire_cut_qpd_frame() -> List[QPDComponent]:
    """
    Constructs the complete 10-term quasiprobability decomposition frame
    for an identity wire cut (gamma = 4.0).

    Expansion identity:
    I(rho) = Tr(Pz+ rho) Pz+ + Tr(Pz- rho) Pz-
           + 0.5 * sum_{s,s' in {+,-}} s*s' * Tr(Px,s rho) Px,s'
           + 0.5 * sum_{s,s' in {+,-}} s*s' * Tr(Py,s rho) Py,s'
    """
    components = [
        # Z-basis direct projections (weight = +1.0)
        QPDComponent(1.0, projector(Z_PLUS), Z_PLUS),
        QPDComponent(1.0, projector(Z_MINUS), Z_MINUS),
        # X-basis cross terms (weight = +/- 0.5)
        QPDComponent(0.5, projector(X_PLUS), X_PLUS),
        QPDComponent(-0.5, projector(X_PLUS), X_MINUS),
        QPDComponent(-0.5, projector(X_MINUS), X_PLUS),
        QPDComponent(0.5, projector(X_MINUS), X_MINUS),
        # Y-basis cross terms (weight = +/- 0.5)
        QPDComponent(0.5, projector(Y_PLUS), Y_PLUS),
        QPDComponent(-0.5, projector(Y_PLUS), Y_MINUS),
        QPDComponent(-0.5, projector(Y_MINUS), Y_PLUS),
        QPDComponent(0.5, projector(Y_MINUS), Y_MINUS),
    ]
    return components


# =====================================================================
# 4. Monolithic Reference Simulator
# =====================================================================

def execute_monolithic_circuit(
    theta1: float, theta2: float, observable: np.ndarray
) -> float:
    """
    Executes the undivided 2-qubit circuit on a single monolithic processor.
    Q0: |0> -> Ry(theta1) -> Rx(theta2) -> Control CNOT -> Measure
    Q1: |0> -----------------------------> Target  CNOT -> Measure
    """
    # Initial state |00>
    psi_00 = np.kron(Z_PLUS, Z_PLUS)

    # Step 1: Ry(theta1) on Q0
    u1 = np.kron(ry_gate(theta1), I2)
    psi_1 = u1 @ psi_00

    # Step 2: Rx(theta2) on Q0
    u2 = np.kron(rx_gate(theta2), I2)
    psi_2 = u2 @ psi_1

    # Step 3: CNOT(Q0 -> Q1)
    psi_out = cnot_gate() @ psi_2

    # Expectation value  = 
    expval = np.real(np.conjugate(psi_out).T @ observable @ psi_out)
    return float(expval)


# =====================================================================
# 5. Distributed QPU Partitions
# =====================================================================

class DistributedQPUCluster:
    """Simulates two isolated QPUs communicating solely via classical channels."""

    @staticmethod
    def run_qpu_a(theta1: float, meas_proj: np.ndarray, shots: int) -> float:
        """
        QPU A: Prepares state, applies Ry(theta1), measures with projector.
        Returns the empirical probability of registering the cut projector.
        """
        psi_init = Z_PLUS
        psi_a = ry_gate(theta1) @ psi_init
        rho_a = projector(psi_a)

        # Exact quantum probability of projector match: Tr(M rho)
        prob = np.real(np.trace(meas_proj @ rho_a))
        prob = np.clip(prob, 0.0, 1.0)

        # Simulate shot measurement statistics
        counts = np.random.binomial(shots, prob)
        return counts / shots

    @staticmethod
    def run_qpu_b(
        prep_state: np.ndarray,
        theta2: float,
        observable: np.ndarray,
        shots: int,
    ) -> float:
        """
        QPU B: Injects prep_state on Q0, initializes Q1 to |0>,
        applies Rx(theta2) on Q0, executes CNOT, and measures observable.
        """
        psi_b_init = np.kron(prep_state, Z_PLUS)
        u_rx = np.kron(rx_gate(theta2), I2)
        psi_b_mid = u_rx @ psi_b_init
        psi_b_final = cnot_gate() @ psi_b_mid

        # Quantum expectation on QPU B
        exact_exp = np.real(
            np.conjugate(psi_b_final).T @ observable @ psi_b_final
        )

        # Add finite-shot measurement variance: Var(O) = 1 - ^2
        variance = max(0.0, 1.0 - exact_exp**2)
        shot_noise = np.random.normal(0.0, np.sqrt(variance / shots))
        return float(np.clip(exact_exp + shot_noise, -1.0, 1.0))


# =====================================================================
# 6. Distributed Execution and Classical Post-Processing
# =====================================================================

def execute_distributed_circuit(
    theta1: float,
    theta2: float,
    observable: np.ndarray,
    shots_per_variant: int,
) -> Tuple[float, float]:
    """
    Coordinates distributed quantum execution via LOCC and QPD reconstruction.
    Returns: (exact_reconstructed_expval, empirical_shot_expval)
    """
    frame = get_wire_cut_qpd_frame()
    exact_reconstructed_val = 0.0
    empirical_sampled_val = 0.0

    for comp in frame:
        # 1. Exact Analytical Subcircuit Pipeline
        psi_a = ry_gate(theta1) @ Z_PLUS
        prob_a_exact = np.real(np.trace(comp.meas_projector @ projector(psi_a)))

        psi_b = cnot_gate() @ (
            np.kron(rx_gate(theta2), I2) @ np.kron(comp.prep_state, Z_PLUS)
        )
        exp_b_exact = np.real(np.conjugate(psi_b).T @ observable @ psi_b)
        exact_reconstructed_val += comp.weight * prob_a_exact * exp_b_exact

        # 2. Distributed Classical LOCC Execution with Finite Shot Noise
        meas_prob_a = DistributedQPUCluster.run_qpu_a(
            theta1, comp.meas_projector, shots_per_variant
        )
        meas_exp_b = DistributedQPUCluster.run_qpu_b(
            comp.prep_state, theta2, observable, shots_per_variant
        )

        # Accumulate weighted classical estimate
        empirical_sampled_val += comp.weight * meas_prob_a * meas_exp_b

    return exact_reconstructed_val, empirical_sampled_val


# =====================================================================
# 7. Verification and Experimentation
# =====================================================================

if __name__ == "__main__":
    np.random.seed(42)

    theta_1 = 0.785398  # pi / 4
    theta_2 = 1.047198  # pi / 3
    shots = 50_000

    observables = [
        ("Z (x) Z", np.kron(PAULI_Z, PAULI_Z)),
        ("X (x) X", np.kron(PAULI_X, PAULI_X)),
        ("Z (x) I", np.kron(PAULI_Z, I2)),
        ("I (x) Z", np.kron(I2, PAULI_Z)),
    ]

    print("=" * 80)
    print("DISTRIBUTED QUANTUM COMPUTING VIA CLASSICAL COMMUNICATION (WIRE CUTTING)")
    print(f"Sampling Overhead per cut (gamma): {sum(abs(c.weight) for c in get_wire_cut_qpd_frame()):.1f}")
    print(f"Shot Multiplier (gamma^2):        {sum(abs(c.weight) for c in get_wire_cut_qpd_frame())**2:.1f}")
    print("=" * 80)
    print(f"{'Observable':<12} | {'Monolithic':<12} | {'QPD Analytical':<15} | {'QPD Distributed (Shots)':<22} | {'Error'}")
    print("-" * 80)

    for name, obs in observables:
        mono_val = execute_monolithic_circuit(theta_1, theta_2, obs)
        exact_qpd, sampled_qpd = execute_distributed_circuit(
            theta_1, theta_2, obs, shots_per_variant=shots
        )
        error = abs(mono_val - sampled_qpd)

        print(
            f"{name:<12} | {mono_val:+11.6f} | {exact_qpd:+14.6f} | "
            f"{sampled_qpd:+21.6f} | {error:.2e}"
        )

    print("=" * 80)
    print("Verification complete: Zero quantum interconnects used.")

5. Hardware Limitations & Future Outlook

While quasiprobability decomposition allows quantum programs to scale across multiple chips today, several physical and algorithmic trade-offs govern its practical deployment.

       100% Coherent Links                                 Pure Classical LOCC
    (Optical Transducers / EPR)                        (Circuit Cutting / QPD)
 ◄─────────────────────────────────────────────────────────────────────────────►
  • Zero shot sampling overhead (γ = 1)               • High shot overhead (γ^(2K))
  • Transducer insertion loss & low fidelity           • Standard classical networking
  • Strict cryogenic optical routing                  • Fully scalable on existing QPUs

The Exponential Sampling Wall

The primary constraint of circuit knitting over classical channels is the exponential growth of the shot multiplier $\gamma_{\text{total}}^2 = \gamma^{2K}$ with respect to the number of cut lines $K$.

For large cut counts ($K > 8$), pure classical decomposition becomes intractable due to finite QPU repetition rates.

Advanced Mitigation Strategies

1. Optimal Frame Design (SIC-POVMs): Replacing standard 6-state Pauli frames with Symmetric Informationally Complete POVMs (e.g., tetrahedral measurement frames) reduces the single-qubit measurement variance and compresses the sampling norm $\gamma$.
2. Maximum-Entanglement Min-Cut Compilers: Circuit graph compilers use integer linear programming (ILP) to place partition boundaries across gates with small rotation angles $\theta \ll \pi/2$, taking advantage of the reduced gate-cut overhead $\gamma(R_{ZZ}(\theta)) = 1 + |\sin(\theta)|$.
3. Hybrid Entanglement-Classical Routing: Near-term architectures will adopt a tiered topology: dense local clusters linked by classical circuit knitting, interconnected across cryostats by sparse, high-fidelity optical EPR links. This hybrid approach limits $K_{\text{classical}}$ while scaling overall qubit counts.

Distributed quantum computing over classical channels transforms inter-chip communication from an intractable hardware challenge into an optimization and classical post-processing workflow. By pairing quasiprobability decomposition with intelligent circuit graph partitioning, modular QPU clusters can execute quantum algorithms well beyond the physical limits of a single cryogenic die.