Noise and decoherence

Density-matrix simulation with noise channels, compared with closed forms, and shot-based estimation as on hardware.

import matplotlib.pyplot as plt
import numpy as np

import metatron_qso as q

A Bell state under depolarising noise

Noise after every gate lowers the purity and the fidelity with the ideal state.

bell = q.Circuit(2).h(0).cx(0, 1)
ideal = q.simulate(bell)
ps = np.linspace(0, 0.3, 31)
fid, pur = [], []
for p in ps:
    rho = q.DensityMatrix.simulate(bell, noise=q.NoiseModel.depolarizing(p, p))
    fid.append(rho.fidelity_with_pure(ideal))
    pur.append(rho.purity())

plt.figure(figsize=(5, 3))
plt.plot(ps, fid, label="fidelity with |Φ⁺⟩")
plt.plot(ps, pur, label="purity")
plt.xlabel("depolarising probability p (one- and two-qubit)")
plt.legend()
plt.tight_layout()
plt.show()

T₁ and T₂

Thermal relaxation of |+⟩ over a time t: the excited population decays as e^{−t/T₁}, the coherence as e^{−t/T₂}.

t1, t2 = 50.0, 30.0
times = np.linspace(0, 150, 61)
pop, coh = [], []
for t in times:
    rho = q.DensityMatrix.simulate(q.Circuit(1).h(0))
    rho = rho.apply_channel(q.Channel.thermal_relaxation(t1, t2, t), [0])
    m = rho.to_matrix()
    pop.append(m[1, 1].real)
    coh.append(2 * abs(m[0, 1]))
print("largest deviation from the closed forms:",
      max(np.max(np.abs(np.array(pop) - 0.5 * np.exp(-times / t1))),
          np.max(np.abs(np.array(coh) - np.exp(-times / t2)))))

plt.figure(figsize=(5, 3))
plt.plot(times, pop, label="excited population")
plt.plot(times, 0.5 * np.exp(-times / t1), "k:", label="½ e^{−t/T₁}")
plt.plot(times, coh, label="coherence 2|ρ₀₁|")
plt.plot(times, np.exp(-times / t2), "k--", label="e^{−t/T₂}")
plt.xlabel("t")
plt.legend()
plt.tight_layout()
plt.show()

largest deviation from the closed forms: 2.220446049250313e-16

Shot noise

estimate_expectation measures groups of commuting Pauli terms with a finite number of shots. The error falls as 1/√shots.

c = q.efficient_su2(4, 2).bind(np.linspace(0.1, 1.5, 24))
h = q.heisenberg_xxz(4, 1.0, 0.5)
exact = h.expectation(q.simulate(c))
shots = [100, 300, 1000, 3000, 10000, 30000]
rms = []
for s in shots:
    est = [q.estimate_expectation(c, h, shots=s, seed=k) for k in range(20)]
    rms.append(np.sqrt(np.mean((np.array(est) - exact) ** 2)))

plt.figure(figsize=(5, 3))
plt.loglog(shots, rms, "o-", label="RMS error (20 seeds)")
plt.loglog(shots, rms[0] * np.sqrt(shots[0] / np.array(shots)), "k:", label="∝ 1/√shots")
plt.xlabel("shots per measurement group")
plt.legend()
plt.tight_layout()
plt.show()

Readout errors

Readout flips act on measured bits only; they do not change the state.

ghz = q.Circuit(3).h(0).cx(0, 1).cx(1, 2)
rho = q.DensityMatrix.simulate(ghz)
print("ideal:      ", rho.sample(4000, seed=1))
print("with readout:", rho.sample(4000, seed=1, readout=(0.02, 0.05)))

ideal:       {'000': 2021, '111': 1979}
with readout: {'000': 1837, '001': 41, '010': 38, '011': 89, '100': 43, '101': 110, '110': 82, '111': 1760}