Dynamic Casimir Photon Emission on a 133-Qubit Quantum Computer

The Centre of mass of excitations above (code on Qwork) shows that the histogram is sharply peaked near index ≈ 32, two to three sites right of the geometric centre (29.5). One σ is about six sites, so ≳ 80% of shots keep their collective excitation weight inside a 1‑D well spanning only 20% of the chain. Dynamical Casimir photons concentrate where the effective mirror crosses the fundamental antinode, but the peak’s right‑shift confirms the asymmetric drive we already saw in the per‑qubit profile. The chip’s control latency (or unequal qubit coupler response) makes the moving boundary settle slightly off‑centre.
The Correlation decay with pair‑pair separation above (code on Qwork) sits above unity everywhere, confirming that symmetric pairs fire more often together than a product model predicts, but it slopes downward from ≈ 1.37. Neighboring mirror‑symmetric pairs share the same standing‑wave lobe, when a photon pair is created at antinode i, the local squeeze field spills into the next few modes, raising their joint‑activation probability by 3 - 5%. Beyond ten pair indices (half a wavelength in the 30-pair basis) the shared lobe overlap drops, so the excess probability erodes toward a baseline. Even the most distant pairs remain ~ 33% more likely to co‑fire than an uncorrelated model (1.0). That plateau reflects a global shot‑to‑shot energy ceiling, shots with high pump energy spur system‑wide pair production, raising all joint counts slightly, while low‑energy shots suppress them. In other words, part of the correlation is collective rather than spatial.
The Pair‑pair covariance above (code on Qwork) shows that positive covariance clusters around pairs 14 - 18 (central antinode) and around 5 - 7 and 23 - 25 (first‑overtone antinodes). Negative covariance appears one step off those hotspots, indicating mode competition, when one pair in a band fires, its neighbors are slightly suppressed, preserving a global energy ceiling shot‑to‑shot. This mirrors photon statistics in a pumped cavity with a fixed total drive energy.
The Energy vs pair multiplicity above (code on Qwork) shows that the 2-D histogram forms a diagonal streak of slope ≈ 11 excited qubits per extra symmetric pair. A single |11〉pair costs ~11 excitations, two pairs ~19 - 20, and so on. Pair creation is accompanied by a nearly constant halo of isolated qubit flips (≈ 5 - 7) that do not participate in symmetric correlations. These could be single‑photon or thermal events generated by the same quench. Efficiency of converting pump energy into correlated photons is about η = 2/11 ≈ 18% (2 qubits out of ~11 carry the Bell‑pair signature).
In the end, using a 60‑qubit chain, this experiment emulated a relativistically moving mirror by applying four sudden ZZ‑quench layers (θ = 0, 0.3π, 0.6π, 0.9π), converted vacuum fluctuations into real excitations, and measured 8 192 hardware shots. The run produced a 72% vacuum‑pair creation rate dominated by single symmetric |11〉pairs, with excitations sharply concentrated near qubits 30‑35, evidence that the effective mirror settled slightly off‑centre due to asymmetric drive latency. Per‑pair statistics revealed mode‑specific multi‑pair squeezing around the central antinode while distance‑averaged correlations remained flat, confirming that entanglement is harmonic, not separation, driven. Energy‑versus‑pair analysis showed ~11 qubit flips per Bell pair, indicating an 18% pump‑to‑correlated‑photon efficiency, and covariance maps exposed competition between neighboring cavity modes under a fixed energy budget. The experiment validates real superconducting qubits as a platform for dynamic Casimir radiation.

Code:

# Main circuit
# Imports 
import json, logging, pandas as pd
from pathlib import Path
from math import pi
from collections import Counter

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit.circuit.library import RZZGate
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2
from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt

logging.basicConfig(level=logging .INFO)
log = logging.getLogger(__name__)

# IBMQ
TOKEN = "IBMQ_KEY_O-`"
INSTANCE = "IBMQ_CRN"
BACKEND_NAME = "ibm_torino"

service = QiskitRuntimeService(
    channel="ibm_cloud",
    token=TOKEN,
    instance=INSTANCE,
)
backend = service.backend(BACKEND_NAME)

# Qubit quality based device mapping 
def best_qubits(csv_path: str, n: int) -> list[int]:
    """Return a list of the `n` best physical qubits by √X‑error → T1 → T2."""
    df = pd .read_csv(csv_path)
    df.columns = df.columns.str.strip()
    order = df.sort_values(
        ["√x (sx) error", "T1 (us)", "T2 (us)"],
        ascending=[True, False, False],
    )
    winners = order["Qubit"].head(n).tolist()
    log .info("Selected %s best qubits: %s", n, winners)
    return winners

# Most recent calibration data
CAL_CSV = Path("/Users/steventippeconnic/Downloads/ibm_torino_calibrations_2025-05-10T17_59_43Z.csv")
PHYSICAL = best_qubits(CAL_CSV, 60)

# Circuit 
q_reg = QuantumRegister(60, "q")
c_reg = ClassicalRegister(60, "c")
qc = QuantumCircuit(q_reg, c_reg, name="DynamicCasimirChain")

# Prepare the vacuum |0…0⟩. Four step 'mirror' schedule (sudden changes in nearest neighbour coupling θ_k).
theta_schedule = [0.0, 0.3 * pi, 0.6 * pi, 0.9 * pi]   
for k, theta in enumerate(theta_schedule):
    log.debug("Applying ZZ layer %d with θ = %.3f rad", k, theta)
    for i in range(59):                                 
        qc.append(RZZGate(theta), [q_reg[i], q_reg[i + 1]])
    qc.barrier(label=f"mirror@step_{k}")

# Measure every qubit
qc.measure(q_reg, c_reg)

# Transpilation 
trans = transpile(
    qc,
    backend=backend,
    initial_layout=PHYSICAL,
    optimization_level=3,
)

#  Execution 
SHOTS = 8192
sampler = SamplerV2(mode=backend)
job = sampler .run([trans], shots=SHOTS)
result = job.result()

creg_name = trans.cregs[0].name  # "c"
counts = result[0].data.__getattribute__(creg_name).get_counts()

#  Derived observable (symmetric‑pair excitations) 
def pair_created(bitstr: str, length: int = 60) -> bool:
    """Return True if *any* symmetric pair (i, L‑1‑i) is in state |11⟩."""
    return any(bitstr[i] == "1" and bitstr[length - 1 - i] == "1" for i in range(length // 2))

pair_events = sum(freq for bits, freq in counts.items() if pair_created(bits))
pair_creation_rate = pair_events / SHOTS
log .info("Vacuum‑pair creation rate: %.5f", pair_creation_rate)

# Correlation count per individual pair 
def pairwise_counts(raw: Counter, length: int = 60):
    pc = Counter()
    for bits, freq in raw.items():
        for i in range(length // 2):
            pair_label = f"{i}-{length-1-i}"
            if bits[i] == bits[length - 1 - i] == "1":
                pc[pair_label] += freq
    return dict(pc)

pairwise = pairwise_counts(counts)

# Json 
out = {
    "experiment_name": "Dynamic Casimir Radiation in Tunable Cavity Chain",
    "theta_schedule_rad": theta_schedule,
    "raw_counts": dict(counts),
    "pair_creation_rate": pair_creation_rate,
    "pairwise_pair_counts": pairwise,
}

JSON_PATH = Path("/Users/steventippeconnic/Documents/Dynamic_Casimir_Cavity_Chain_0.json")
with JSON_PATH.open("w") as fp:
    json.dump(out, fp, indent=4)
log .info("Results saved → %s", JSON_PATH)

#  Visual
plot_histogram(counts, title="Dynamic Casimir Chain – Measured Bitstrings")
plt .show()

# End







/////////////////////////////////////////////////////////////////
# Code for all visuals from experiment JSON

# Imports
import json
from pathlib import Path
from collections import Counter
import numpy as np
import matplotlib.pyplot as plt

# Load data
file_path = Path("/Users/steventippeconnic/Documents/Dynamic_Casimir_Cavity_Chain_0.json")
data = json.loads(file_path.read_text())
counts_raw = Counter(data["raw_counts"])
chain_len = 60
shots = sum(counts_raw.values())
pairs = [(i, chain_len - 1 - i) for i in range(chain_len // 2)]
L = 60

# Prepare histograms and tallies
pairs_per_shot = Counter()
hamming_hist = Counter()
pair_counts = Counter()
pair_corr = np.zeros((chain_len // 2, chain_len // 2))

for bits, freq in counts_raw.items():
    active_pairs = []
    for idx, (i, j) in enumerate(pairs):
        if bits[i] == bits[j] == "1":
            pair_counts[idx] += freq
            active_pairs.append(idx)
    pairs_per_shot[len(active_pairs)] += freq
    for a in active_pairs:
        pair_corr[a, a] += freq
        for b in active_pairs:
            if a != b:
                pair_corr[a, b] += freq
    hamming_hist[bits.count("1")] += freq

pair_corr /= shots

# Pair multiplicity histogram
plt.figure()
plt.bar(pairs_per_shot.keys(), pairs_per_shot.values())
plt.xlabel("Symmetric pair count per shot")
plt.ylabel("Frequency")
plt.title("Pair multiplicity distribution")
plt.show()

# Spatial pair distribution
plt.figure()
indices, values = zip(*sorted(pair_counts.items()))
plt.bar(indices, values)
plt.xlabel("Pair index i and L minus 1 minus i")
plt.ylabel("Total events")
plt.title("Spatial distribution of created pairs")
plt.show()

# Pair correlation heat map
plt.figure()
plt.imshow(pair_corr, aspect="auto", interpolation="nearest")
plt.colorbar(label="Joint probability")
plt.xlabel("Pair j")
plt.ylabel("Pair i")
plt.title("Pair to pair correlation matrix")
plt.show()

# Hamming weight histogram
plt.figure()
plt.bar(hamming_hist.keys(), hamming_hist.values(), width=0.8)
plt.xlabel("Total excited qubits in shot")
plt.ylabel("Frequency")
plt.title("Excitation energy distribution")
plt.show()

# Helper structures
excitation_prob = np.zeros(L)
longest_cluster_hist = Counter()
side_diff_hist = Counter()
pair_joint = np.zeros((L // 2, L // 2))
pair_single = np.zeros(L // 2)
pairs = [(i, L - 1 - i) for i in range(L // 2)]

# Analyse each shot
for bitstring, freq in counts_raw.items():
    bits = np.array(list(bitstring), dtype=int)
    excitation_prob += freq * bits

    # Longest consecutive cluster of ones
    runs = np.diff(np.where(np.concatenate(([bits[0]],
                                            bits[:-1] ^ bits[1:],
                                            [1])))[0])[::2]
    longest_cluster_hist[max(runs)] += freq

    # Left‑minus‑right excitation difference
    left_ones = bits[:L // 2].sum()
    right_ones = bits[L // 2:].sum()
    side_diff_hist[left_ones - right_ones] += freq

    # Pair statistics
    active_pairs = []
    for idx, (i, j) in enumerate(pairs):
        if bits[i] == bits[j] == 1:
            pair_single[idx] += freq
            active_pairs.append(idx)
    for i in active_pairs:
        for j in active_pairs:
            pair_joint[i, j] += freq

excitation_prob /= shots
pair_single /= shots
pair_conditional = np.divide(pair_joint / shots,
                             pair_single[:, None],
                             out=np.zeros_like(pair_joint),
                             where=pair_single[:, None] > 0)

# Per‑qubit excitation probability
plt.figure()
plt.plot(range(L), excitation_prob, marker="o")
plt.xlabel("Qubit index")
plt.ylabel("Probability qubit is 1")
plt.title("Per‑qubit excitation profile")
plt.show()

# Histogram of longest 1‑run per shot
plt.figure()
plt.bar(longest_cluster_hist.keys(), longest_cluster_hist.values())
plt.xlabel("Longest consecutive 1‑run length")
plt.ylabel("Frequency")
plt.title("Cluster size distribution")
plt.show()

# Conditional pair activation heat map
plt.figure()
plt.imshow(pair_conditional, aspect="auto", interpolation="nearest")
plt.colorbar(label="P(pair j | pair i)")
plt.xlabel("Pair j")
plt.ylabel("Pair i")
plt.title("Conditional activation probability matrix")
plt.show()

# Left‑minus‑right excitation balance
plt.figure()
plt.bar(side_diff_hist.keys(), side_diff_hist.values(), width=0.8)
plt.xlabel("Left half minus right half excitations")
plt.ylabel("Frequency")
plt.title("Spatial asymmetry of excitations")
plt.show()

# Accumulators
pair_single = np.zeros(L // 2)
pair_joint = np.zeros((L // 2, L // 2))
dist_accum = np.zeros(L // 2)
dist_counts = np.zeros(L // 2)
com_hist = Counter()
pair_energy_hist = Counter()

# Main sweep
for bits, freq in counts_raw.items():
    bit_arr = np.array([ch == "1" for ch in bits], dtype=int)
    active_pairs = [p_idx
                    for p_idx, (i, j) in enumerate(pairs)
                    if bit_arr[i] and bit_arr[j]]
    
    #   Joint accumulator one increment per ordered pair
    for a in active_pairs:
        for b in active_pairs:
            d = abs(a - b)
            dist_accum[d] += freq            

    #   Denominator one increment per *distance* in this shot
    seen_d = {abs(a - b) for a in active_pairs for b in active_pairs}
    for d in seen_d:
        dist_counts[d] += freq       

    # Centre of mass
    idx_excited = np.flatnonzero(bit_arr)
    if idx_excited.size:
        com_hist[int(round(idx_excited.mean()))] += freq

    # Symmetric pair data
    active_pairs = []
    for p_idx, (i, j) in enumerate(pairs):
        if bit_arr[i] and bit_arr[j]:
            pair_single[p_idx] += freq
            active_pairs.append(p_idx)

    # Energy vs pair count
    pair_energy_hist[(len(active_pairs), bit_arr.sum())] += freq

    # Correlations by separation
    for a in active_pairs:
        for b in active_pairs:
            pair_joint[a, b] += freq
            d = abs(a - b)
            dist_accum[d] += freq
            dist_counts[d] += freq           

pair_single /= shots
pair_joint /= shots
cov = pair_joint - np.outer(pair_single, pair_single)

# Centre‑of‑mass histogram
plt.figure()
if com_hist:
    x, y = zip(*sorted(com_hist.items()))
    plt.bar(x, y)
plt.xlabel("Excitation centre of mass")
plt.ylabel("Frequency")
plt.title("Centre of mass of excitations")
plt.show()

# Correlation vs separation
sep = np.arange(1, L // 2)
corr_vs_d = np.where(dist_counts[1:] > 0,
                     dist_accum[1:] / dist_counts[1:],
                     0.0)
plt.figure()
plt.plot(sep, corr_vs_d, marker="o")
plt.xlabel("Pair‑pair separation")
plt.ylabel("Average joint probability")
plt.title("Correlation decay with separation")
plt.show()

# Covariance heat map
plt.figure()
plt.imshow(cov, cmap="seismic", aspect="auto", interpolation="nearest")
plt.colorbar(label="Covariance")
plt.xlabel("Pair j")
plt.ylabel("Pair i")
plt.title("Pair‑pair covariance")
plt.show()

# 2‑D histogram pair count vs energy
max_pairs = max(pc for pc, _ in pair_energy_hist.keys())
max_energy = max(en for _, en in pair_energy_hist.keys())
grid = np.zeros((max_pairs + 1, max_energy + 1))
for (pc, en), f in pair_energy_hist.items():
    grid[pc, en] = f
plt.figure()
plt.imshow(grid, origin="lower", aspect="auto", interpolation="nearest")
plt.colorbar(label="Frequency")
plt.xlabel("Total excited qubits")
plt.ylabel("Number of symmetric pairs")
plt.title("Energy vs pair multiplicity")
plt.show()

# End