We're trying to simulate erasure errors on the surface code using Stim. The threshold for erasure errors on the data qubits (after initialization) is 50%.
We followed the following post: How do I perform an erasure error in stim?
However, we end up with a threshold of ~20%. This is the same as this post: Threshold value when simulating erasures with stim
That post did not get an answer and so we bring it up again here.
This is our circuit for the unrotated $d=3$ surface code
RX 0 2 4 10 12 14 20 22 24 6 8 16 18 1 3 5 7 9 11 13 15 17 19 21 23
TICK
R 999
E(0.025) X999
E(0.025) X0 X999
E(0.025) Y0 X999
E(0.025) Z0 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X2 X999
E(0.025) Y2 X999
E(0.025) Z2 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X4 X999
E(0.025) Y4 X999
E(0.025) Z4 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X10 X999
E(0.025) Y10 X999
E(0.025) Z10 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X12 X999
E(0.025) Y12 X999
E(0.025) Z12 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X14 X999
E(0.025) Y14 X999
E(0.025) Z14 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X20 X999
E(0.025) Y20 X999
E(0.025) Z20 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X22 X999
E(0.025) Y22 X999
E(0.025) Z22 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X24 X999
E(0.025) Y24 X999
E(0.025) Z24 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X6 X999
E(0.025) Y6 X999
E(0.025) Z6 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X8 X999
E(0.025) Y8 X999
E(0.025) Z8 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X16 X999
E(0.025) Y16 X999
E(0.025) Z16 X999
M 999
DETECTOR rec[-1]
R 999
E(0.025) X999
E(0.025) X18 X999
E(0.025) Y18 X999
E(0.025) Z18 X999
M 999
DETECTOR rec[-1]
CZ 6 1 8 3 16 11 18 13
CX 5 10 7 12 9 14 15 20 17 22 19 24
TICK
CZ 2 1 4 3 12 11 14 13 22 21 24 23
CX 5 6 7 8 15 16 17 18
TICK
CZ 0 1 2 3 10 11 12 13 20 21 22 23
CX 7 6 9 8 17 16 19 18
TICK
CZ 6 11 8 13 16 21 18 23
CX 5 0 7 2 9 4 15 10 17 12 19 14
TICK
MRX 1 3 5 7 9 11 13 15 17 19 21 23
DETECTOR rec[-10]
DETECTOR rec[-9]
DETECTOR rec[-8]
DETECTOR rec[-5]
DETECTOR rec[-4]
DETECTOR rec[-3]
REPEAT 2 {
TICK
CZ 6 1 8 3 16 11 18 13
CX 5 10 7 12 9 14 15 20 17 22 19 24
TICK
CZ 2 1 4 3 12 11 14 13 22 21 24 23
CX 5 6 7 8 15 16 17 18
TICK
CZ 0 1 2 3 10 11 12 13 20 21 22 23
CX 7 6 9 8 17 16 19 18
TICK
CZ 6 11 8 13 16 21 18 23
CX 5 0 7 2 9 4 15 10 17 12 19 14
TICK
MRX 1 3 5 7 9 11 13 15 17 19 21 23
DETECTOR rec[-1] rec[-53]
DETECTOR rec[-2] rec[-54]
DETECTOR rec[-3] rec[-55]
DETECTOR rec[-4] rec[-56]
DETECTOR rec[-5] rec[-57]
DETECTOR rec[-6] rec[-58]
DETECTOR rec[-7] rec[-59]
DETECTOR rec[-8] rec[-60]
DETECTOR rec[-9] rec[-61]
DETECTOR rec[-10] rec[-62]
DETECTOR rec[-11] rec[-63]
DETECTOR rec[-12] rec[-64]
}
MX 0 2 4 6 8 10 12 14 16 18 20 22 24
DETECTOR rec[-23] rec[-8] rec[-10] rec[-13]
DETECTOR rec[-22] rec[-7] rec[-9] rec[-10] rec[-12]
DETECTOR rec[-21] rec[-6] rec[-9] rec[-11]
DETECTOR rec[-18] rec[-3] rec[-5] rec[-8]
DETECTOR rec[-17] rec[-2] rec[-4] rec[-5] rec[-7]
DETECTOR rec[-16] rec[-1] rec[-4] rec[-6]
OBSERVABLE_INCLUDE(0) rec[-13] rec[-12] rec[-11]
Note: If we simulate our codes under pauli data noise we recover the expected thresholds so the error is not in the circuit for the codes but instead in the way we are handling erasures.
Decoding
Decoding is done using Sinter to call PyMatching, as shown in the Stim tutorial. We passed 2 new flags to the detector error model: approximate_disjoint_errors=True and ignore_decomposition_failures=True. See the code below.
import stim
import numpy as np
import pymatching
import sinter
def count_logical_errors(circuit: stim.Circuit, num_shots: int) -> int:
# Sample the circuit.
sampler = circuit.compile_detector_sampler()
detection_events, observable_flips = sampler.sample(num_shots, separate_observables=True)
# Extract decoder configuration data from the circuit.
detector_error_model = circuit.detector_error_model(decompose_errors=True,
approximate_disjoint_errors=True,
ignore_decomposition_failures=True)
# Run the decoder.
predictions = sinter.predict_observables(
dem=detector_error_model,
dets=detection_events,
decoder='pymatching',
)
# Count the mistakes.
num_errors = 0
for actual_flip, predicted_flip in zip(observable_flips, predictions):
if not np.array_equal(actual_flip, predicted_flip):
num_errors += 1
return num_errors
```
