Density-matrix simulation of small surface codes under current and projected experimental noise

We present a density-matrix simulation of the quantum memory and computing performance of the distance-3 logical qubit Surface-17, following a recently proposed quantum circuit and using experimental error parameters for transmon qubits in a planar circuit QED architecture. We use this simulation to optimize components of the QEC scheme (e.g., trading off stabilizer measurement infidelity for reduced cycle time) and to investigate the benefits of feedback harnessing the fundamental asymmetry of relaxation-dominated error in the constituent transmons. A lower-order approximate calculation extends these predictions to the distance-5 Surface-49. These results clearly indicate error rates below the fault-tolerance threshold of the surface code, and the potential for Surface-17 to perform beyond the break-even point of quantum memory. However, Surface-49 is required to surpass the break-even point of computation at state-of-the-art qubit relaxation times and readout speeds. Quantum computing hardware has reached low enough error rates for any computation to be performed given a large enough system. A collaboration between Delft University of Technology and Leiden University have simulated the performance of state-of-the-art superconducting qubits in a small quantum circuit designed to protect information from otherwise catastrophic errors. Such a circuit requires the underlying device to perform above a ‘fault-tolerant’ threshold before it can be scaled to a full quantum computer. The team determined that not only should current technology break this barrier, but also that this quantum circuit should already store information better than its individual components. This work provides performance targets and guidance on design choices for experiments currently in development across the world.


INTRODUCTION
Recent experimental demonstrations of small quantum simulations 1-3 and quantum error correction (QEC) 4-7 position superconducting circuits for targeting quantum supremacy 8 and quantum fault tolerance, 9 two outstanding challenges for all quantum information processing platforms. On the theoretical side, much modeling of QEC codes has been made to determine fault-tolerance threshold rates in various models [10][11][12] with different error decoders. [13][14][15] However, the need for computational efficiency has constrained many previous studies to oversimplified noise models, such as depolarizing and bit-flip noise channels. This discrepancy between theoretical descriptions and experimental reality compromises the ability to predict the performance of near-term QEC implementations, and offers limited guidance to the experimentalist through the maze of parameter choices and trade-offs. In the planar circuit quantum electrodynamics (cQED) 16 architecture, the major contributions to error are transmon qubit relaxation, dephasing from flux noise, and resonator photons leftover from measurement, and leakage from the computational space, none of which are wellapproximated by depolarizing or bit-flip channels. Simulations with more complex error models are now essential to accurately pinpoint the leading contributions to the logical error rate in the small-distance surface codes 10,13,17 currently pursued by several groups worldwide.
In this paper, we perform a density-matrix simulation of the distance-3 surface code named Surface-17, using the concrete quantum circuit recently proposed in ref. 18 and the measured performance of current experimental multi-transmon cQED platforms. [19][20][21][22] For this purpose, we have developed an open-source density-matrix simulation package named quantumsim (https:// github.com/brianzi/quantumsim). We use quantumsim to extract the logical error rate per QEC cycle, ϵ L . This metric allows us to optimize and trade-off between QEC cycle parameters, assess the merits of feedback control, predict gains from future improvements in physical qubit performance, and quantify decoder performance. We compare an algorithmic decoder using minimum-weight perfect matching (MWPM) with homemade weight calculation to a simple look-up table (LT) decoder, and weigh both against an upper bound (UB) for decoder performance obtainable from the density-matrix simulation. Finally, we make a low-order approximation to extend our predictions to the distance-5 Surface-49. The combination of results for Surface-17 and Surface-49 allows us to make statements about code scaling and to predict the code size and physical qubit performance required to achieve break-even points for memory and computational performance.

RESULTS
Error rates for Surface-17 under current experimental conditions To quantify the performance of the logical qubit, we first define a test experiment to simulate. Inspired by the recent experimental demonstration of distance-3 and distance-5 repetition codes, 4 we first focus on the performance of the logical qubit as a quantum memory. Specifically, we quantify the ability to hold a logical 0 j i state, by initializing this state, holding it for k ∈ {1, …, 20} cycles, performing error correction, and determining a final logical state (see Fig. 6 for details). The logical fidelity F L [k] is then given by the probability to match the initial state. We observe identical results when using 1 j i or ± j i ¼ 1 ffiffi Þin place of 0 j i. We base our error model for the physical qubits on current typical experimental performance for transmons in planar cQED, using parameters from the literature and in-house results (e.g., gate-set tomography measurements). These are summarized in Table 1, and further detailed in Supplemental Material. We focus on the QEC cycle proposed in ref. 18 which pipelines the execution of X-type and Z-type stabilizer measurements. Each stabilizer measurement consists of three parts: a coherent step (duration τ c = 2τ g,1Q + 4τ g,2Q ), measurement (τ m ), and photon depletion from readout resonators (τ d ), making the QEC cycle time τ cycle = τ c + τ m + τ d .
Simulating this concrete quantum circuit with the listed parameters using quantumsim, we predict F L [k] of Surface-17 ( Fig. 1). We show F L [k] for both a homemade MWPM decoder (green, described in Supplemental Material), and an implementation of the LT decoder of ref. 13 (blue, described in Supplemental Material). To isolate decoder performance, we can compare the achieved fidelity to an UB extractable from the density-matrix simulation (red, described in Methods). To assess the benefit of QEC, we also compare to a single decohering transmon, whose fidelity is calculated by averaging over the six cardinal points of the Bloch sphere: The observation of F L [k] > F phys (kτ cycle ) for large k would constitute a demonstration of QEC beyond the quantum memory break-even point. 7 Equivalently, one can extract a logical error rate ϵ L from a best fit to F L [k] (as derived in Methods as the probability of an odd number of errors occurring), Here, k 0 and ϵ L are the parameters to be fit. We compare ϵ L to the physical error rate We observe ϵ L = 1.44% c for the LT decoder, ϵ L = 1.07% c for the MWPM decoder, and ϵ L = 0.68% c at the decoder UB (% c = % per cycle). The latter two fall below ϵ phys = 1.33% c . Defining the decoder efficiency We can also compare the multi-cycle error correction to majority voting, in which the state declaration is based solely on the output of the final data qubit measurements (ancilla measurements are ignored). Majority voting corrects any single data qubit error (over the entire experiment), and thus exhibits a quadratic decay for small k (A distance-d code with majority voting alone should exhibit a (d + 1)/2-order decay). A decoder should also be able to correct (at least) a single error, and thus should produce the same behavior at low k, delaying the onset of exponential decay in F L [k]. In fact, a good test for the performance of a MWPM decoder is to ensure it can outperform the majority vote at short timescales, as suboptimal configuration will prevent this (as seen for the LT decoder).
With the baseline for current performance established, we next investigate ϵ L improvements that may be achieved by two means. First, we consider modifications to the QEC cycle at fixed physical performance. Afterwards, we consider the effect of improving physical qubit T 1 and T ϕ .
Optimization of logical error rates with current experimental conditions Error sources in current cQED setups derive primarily from transmon decoherence, as opposed to gate and measurement errors produced by control electronics. Thus, a path to reducing ϵ L may be to decrease τ cycle . Currently, the cycle is dominated by τ m + τ d . At fixed readout power, reducing τ m and τ d will reduce τ cycle at the cost of increased readout infidelity ϵ RO (described in Methods). We explore this trade-off in Fig. 2, using a lineardispersive readout model, 23 keeping τ m = τ d and assuming no leftover photons. Because of the latter, ϵ  is smaller for τ m = 200 ns than for τ m = 300 ns, even though ϵ RO increases to 5%. It is interesting to note that the optimal τ m for quantum memory, which minimizes logical error per unit time, rather than per cycle, is τ m = 280 ns (Fig. 2 inset). This shows that different cycle parameters might be optimal for computation and memory applications.
Next, we consider the possibility to reduce ϵ L using feedback control. Since T 1 only affects qubits in the excited state, the error rate of ancillas in Surface-17 is roughly two times higher when in the excited state. The unmodified syndrome extraction circuit flips the ancilla if the corresponding stabilizer value is −1, and since ancillas are not reset between cycles, they will spend significant amounts of time in the excited state. Thus, we consider using feedback to hold each ancilla in the ground state as much as possible. We do not consider feedback on data qubits, as the highly entangled logical states are equally susceptible to T 1 .
The feedback scheme (Inset of Fig. 3) consists of replacing the R y (π/2) gate at the end of the coherent step with a R y (−π/2) gate for some of the ancillas, depending on a classical control bit p for each ancilla. This bit p represents an estimate of the stabilizer value, and the ancilla is held in the ground state whenever this estimate is correct (i.e., in the absence of errors). Figure 3 shows the effect of this feedback on the logical fidelity, both for the MWPM decoder and the decoder UB. We observe ϵ L improve only 0.05% c in both cases. Future experiments might opt not to pursue these small gains in view of the technical challenges added by feedback control.
Projected improvement with advances in quantum hardware We now estimate the performance increase that may result from improving the transmon relaxation and dephasing times via materials and filtering improvements. To model this, we return to τ cycle = 800 ns, and adjust T 1 values with both T ϕ = 2T 1 (common in experiment) and T ϕ = ∞ (all white-noise dephasing eliminated). We retain the same rates for coherent errors, readout infidelity, and photon-induced dephasing as in Fig. 1. Figure 4 shows the extracted ϵ L and ϵ phys over the T 1 range covered. For the MWPM decoder (UB) and T ϕ = 2T 1 , the memory figure of merit γ m = ϵ phys /ϵ L increases from 1.3 (2) at T 1 = 30 μs to 2 (5) at 100 μs.
Completely eliminating white-noise dephasing will increase γ m by 10% with MWPM and 30% at the UB.
A key question for any QEC code is how ϵ L scales with code distance d. Computing power limitations preclude similar densitymatrix simulations of the d = 5 surface code Surface-49. However, we can approximate the error rate by summing up all lowest-order error chains (as calculated for the MWPM decoder), and deciding individually whether or not these would be corrected by a MWPM decoder (See Supplemental Material for details). Figure 5 shows the lowest-order approximation to the logical error rates of Surface-17 and Surface-49 over a range of T 1 = T ϕ /2. Comparing the Surface-17 lowest-order approximation to the quantumsim result shows good agreement and validates the approximation. We observe a lower ϵ L for Surface-49 than for Surface-17, indicating quantum fault tolerance over the T 1 range covered. The fault-tolerance figure of merit defined in, 9 increases from 2 to 4 as T 1 grows from 30 to 100 μs.
As a rough metric of computational performance, we offer to compare ϵ L (per cycle) to the error accrued by a physical qubit Fig. 2 Optimization of the logical error rate (per cycle) of Surface-17 as a function of measurement-and-depletion time. 19 Changes in the underlying physical error rates are shown, as well. Decreasing the measurement time causes an increase in the readout infidelity (solid black curve with dots), whilst decreasing the single qubit decay from T 1 and T 2 (black dashed curve) for all qubits. The logical rate with an MWPM decoder (green curve) is minimized when these error rates are appropriately balanced. The logical error rate is calculated from the best fit of Eq. (2). Error bars (2 s.d.) are obtained by bootstrapping (N = 10,000 runs). Inset: Logical error rate per unit time, instead of per cycle  2) to the data, and error bars (2 s.d.) are given by bootstrapping, with each point averaged over 10,000 runs. Inset: Method for implementing the feedback scheme. For each ancilla qubit A j , we store a parity bit p j , which decides the sign of the R y (π/2) rotation at the end of each coherent step. The time A j spends in the ground state is maximized when p j is updated each cycle t by XORing with the measurement result from cycle t − 1, after the rotation of cycle t has been performed idling over τ g,1Q . We define a metric for computation performance, and γ c = 1 as a computational breakeven point. Clearly, using the QEC cycle parameters of Table 1 and even with T 1 improvements, neither Surface-17 nor Surface-49 can break-even computationally. However, including the readout acceleration recently demonstrated in ref. 22 which allows τ m = τ d = 100 ns and τ cycle = 400 ns, Surface-49 can cross γ c = 1 by T 1 = 40 μs . In view of first reports of T 1 up to 80 μs emerging for planar transmons, 24,25 this important milestone may be within grasp.

DISCUSSION
Computational figure of merit We note that our metric of computational power is not rigorous, due to the different gate sets available to physical and logical qubits. Logical qubits can execute multiple logical X and Z gates within one QEC cycle, but require a few cycles for two-qubit and Hadamard gates (using the proposals of refs. 12, 17), and state distillation over many cycles to perform non-Clifford gates. As such, this metric is merely a rough benchmark for computational competitiveness of the QEC code. However, given the amount by which all distance-3 logical fidelities fall above this metric, we find it unlikely that these codes will outperform a physical qubit by any fair comparison in the near future.
Decoder performance A practical question facing QEC is how best to balance the tradeoff between decoder complexity and performance. Past proposals for surface-code computation via lattice surgery 17 require the decoder to provide an up-to-date estimate of the Pauli error on physical qubits during each logical T gate. Because tracking Pauli errors through a non-Clifford gate is inefficient, however implemented, equivalent requirements will hold for any QEC code. 26 A decoder is thus required to process ancilla measurements from one cycle within the next (on average). This presents a considerable challenge for transmon-cQED implementations, as τ cycle < 1 μs. This short time makes the use of computationally intensive decoding schemes difficult, even if they provide lower ϵ L The leading strategy for decoding the surface code is MWPM using the blossom algorithm of Edmonds. 10,14,27 Although this algorithm is challenging to implement, it scales linearly in code distance. 27 The algorithm requires a set of weights (representing the probability that two given error signals are connected by a chain of errors) as input. An important practical question (See Supplemental Material) is whether these weights can be calculated on the fly, or must be precalculated and stored. On-thefly weight calculation is more flexible. For example, it can take into account the difference in error rates between an ancilla measured in the ground and in the excited state. The main weakness of MWPM is the inability to explicitly detect Y errors. In fact, in Supplemental Material we show that MWPM is nearly perfect in the absence of Y errors. The decoder efficiency Y may significantly increase by extending MWPM to account for correlations between detected X and Z errors originating from Y errors. 28,29 If computational limitations preclude a MWPM decoder from keeping up with τ cycle , the LT decoder may provide a straightforward solution for Surface-17. However, at current physical performance, the η d reduction will make Surface-17 barely miss memory break-even (Fig. 1). Furthermore, memory requirements make LT decoding already impractical for Surface-49. Evidently, real-time algorithmic decoding by MWPM or improved variants is an important research direction already at low code distance.

Other observations
The simulation results allow some further observations. Although we have focused on superconducting qubits, we surmise that the following statements are fairly general.
We observe that small quasi-static qubit errors are suppressed by the repeated measurement. In our simulations, the 1/f flux noise producing 0.01 radians of phase error per flux pulse on a qubit has a diamond norm approximately equal to T 1 noise, but a trace distance 100 times smaller. As the flux noise increases ϵ L by only 0.01% c , it appears ϵ L is dependent on the trace distance rather than the diamond norm of the underlying noise components. Quasi-static qubit errors can then be easily suppressed, but will also easily poison an experiment if unchecked.
We further observe that above a certain value, ancilla and measurement errors have a diminished effect on ϵ L . In our error model, the leading sources of error for a distance d code are chains of (d − 1)/2 data qubit errors plus either a single ancilla qubit error or readout error, which together present the same syndrome as a chain of (d + 1)/2 data qubit errors. An optimal decoder decides which of these chains is more likely, at which point the less-likely chain will be wrongly corrected, completing a logical error. This implies that if readout infidelity ϵ RO ð Þ or the ancilla error rate ϵ anc ð Þ is below the data qubit ϵ phys À Á error rate, . However, if ϵ RO ϵ anc ð Þ>ϵ phys , ϵ L becomes independent of ϵ RO ϵ anc ð Þ, to lowest order. This can be seen in Fig. 2, where the error rate is almost constant as ϵ RO exponentially increases. This approximation breaks down with large enough ϵ anc and ϵ RO , but presents a counter-intuitive point for experimental design; ϵ L becomes less sensitive to measurement and ancilla errors as these error get worse.
A final, interesting point for future surface-code computation is shown in Fig. 2: the optimal cycle parameters for logical error rates per cycle and per unit time are not the same. This implies that logical qubits functioning as a quantum memory should be treated differently to those being used for computation. This idea can be extended further: at any point in time, a large quantum computer performing a computation will have a set S m of memory qubits which are storing part of a large entangled state, whilst a set S c of computation qubits containing the rest of the state undergo operations. To minimize the probability of a logical error occurring on qubits within both S c and S m , the cycle time of the qubits in S c can be reduced to minimize the rest time of qubits in S m . As a simple example, consider a single computational qubit q c and a single memory qubit q m sharing entanglement. Operating all qubits at τ cycle = 720 ns to minimize ϵ L would lead to a 1.09% error rate for the two qubits combined. However, shortening the τ cycle of q c reduces the time over which q m decays. If q c operates at τ cycle = 600 ns, the average error per computational cycle drops to 1.06%, as q m completes only 5 cycles for every 6 on q c . Although this is only a meager improvement, one can imagine that when many more qubits are resting than performing computation, the relative gain will be quite significant.
Effects not taken into account Although we have attempted to be thorough in the detailing of the circuit, we have neglected certain effects. We have used a simple model for C-Z gate errors as we lack data from experimental tomography (e.g., one obtained from two-qubit gate-set tomography 30 ). Most importantly, we have neglected leakage, where a transmon is excited out of the two lowest energy states, i.e., out of the computational subspace. Previous experiments have reduced the leakage probability per C-Z gate tõ 0.3%, 31 and per single-qubit gate to~0.001%. 32 Schemes have also been developed to reduce the accumulation of leakage. 33 Extending quantumsim to include and investigate leakage is a next target. However, the representation of the additional quantum state can increase the simulation effort significantly [by a factor of (9/4) 10 ≈ 3000]. To still achieve this goal, some further approximations or modifications to the simulation will be necessary. Future simulations will also investigate the effect of spread in qubit parameters, both in space (i.e., variation of physical error rates between qubits) and time (e.g., T 1 fluctuations), and cross-talk effects such as residual couplings between nearest and next-nearest neighbor transmons, qubit cross-driving, and qubit dephasing by measurement pulses targeting other qubits.

Simulated experimental procedure
Surface-17 basics. A QEC code can be defined by listing the data qubits and the stabilizer measurements that are repeatedly performed upon them. 34 In this way, Surface-17 is defined by a 3 × 3 grid of data qubits {D 0 , … D 8 }. In order to stabilize a single logical qubit, 9 − 1 = 8 commuting measurements are performed. The stabilizers are the weight-two and weight-four X-type and Z-type parity operators X 2 X 1 , , and X 7 X 6 , where X j (Z j ) denotes the X (Z) Pauli operator acting on data qubit D j . Their measurement is realized indirectly using nearest-neighbor interactions between data and ancilla qubits arranged in a square lattice, followed by ancilla measurements [Fig. 6a]. This leads to a total of 17 physical qubits when a separate ancilla is used for each individual measurement. We follow the circuit realization of this code described in ref. 18, for which we give a schematic description in Fig. 6b (See Supplemental Material for a full circuit diagram).
In an experimental realization of this circuit, qubits will regularly accumulate errors. Multiple errors that occur within a short period of time Fig. 6 Schematic overview of the simulated experiment. a 17 qubits are arranged in a surface code layout (legend top-right). The red data qubits are initialized in the ground state 0 j i, and projected into an eigenstate of the measured X-(blue) and Z-(green) type stabilizer operators. b A section of the quantum circuit depicting the four-bit parity measurement implemented by the A 3 ancilla qubit (+/− refer to R y (±π/2) single-qubit rotations). The ancilla qubit (green line, middle) is entangled with the four data qubits (red lines) to measure Z 1 Z 2 Z 4 Z 5 . Ancillas are not reset between cycles. Instead, the implementation relies on the quantum non-demolition nature of measurements. The stabilizer is then the product of the ancilla measurement results of successive cycles. This circuit is performed for all ancillas and repeated k times before a final measurement of all (data and ancilla) qubits. c All syndrome measurements of the k cycles are processed by the decoder. d After each cycle, the decoder updates its internal state to represent the most likely set of errors that occurred. e After the final measurement, the decoder uses the readout from the data qubits, along with previous syndrome measurements, to declare a final logical state. To this end, the decoder processes the Z-stabilizers obtained directly from the data qubits, finalizing its prediction of most likely errors. The logical parity is then determined as the product of all data qubit parities Q 8 j¼0 D j once the declared errors are corrected. The logical fidelity F L is the probability that this declaration is the same as the initial state 0 j i ð Þ Density-matrix simulation of small surface TE O'Brien et al.
(e.g., one cycle) form error 'chains' that spread across the surface. Errors on single qubits, or correlated errors within a small subregion of Surface-17, fail to commute with the stabilizer measurements, creating error signals that allow diagnosis and correction of the error via a decoder. However, errors that spread across more than half the surface in a short enough period of time are misdiagnosed, causing an error on the logical qubit when wrongly corrected. 10 The rate at which these logical errors arise is the main focus of this paper.
Protocol for measurement of logical error rates. As the performance measure of Surface-17, we study the fidelity of the logical qubit as a quantum memory. We describe our protocol with an example 'run' in Fig. 6. We initialize all qubits in 0 j i and perform k = 1, 2, …, 20 QEC cycles [ Fig. 6b]. Although this initial state is not a stabilizer eigenstate, the first QEC cycle projects the system into one of the 16 overlapping eigenstates within the +1 eigenspace for Z stabilizers, which form the logical 0 j i state. 10 This implies that, in the absence of errors, the first measurement of the Z stabilizers will be +1, whilst that of the X stabilizers will be random. In the following cycles, ancilla measurements of each run [ Fig. 6c] are processed using a classical decoding algorithm. The decoder computes a Pauli update after each QEC cycle [ Fig. 6d]. This is a best estimate of the Pauli operators that must be applied to the data qubits to transform the logical qubit back to the logical 0 j i state. The run ends with a final measurement of all data qubits in the computational basis. From this 9-bit outcome, a logical measurement result is declared [ Fig. 6e]. First, the four Z-type parities are calculated from the 9 data-qubit measurement outcomes and presented to the decoder as a final set of parity measurements. This ensures that the final computed Pauli update will transform the measurement results into a set that measures +1 for all Z stabilizers. This results in one of 32 final measurements, from which the value of a logical Z operator can be calculated to give the measurement result (any choice of logical operator gives the same result). The logical fidelity F L [k] after k QEC cycles is defined as the probability of this declared result matching the initial +1 state.
At long times and with low error rates, Surface codes have a constant logical error rate ϵ L . The fidelity F L [k] is obtained by counting the probability of an odd number of errors having occurred in total (as two σ x errors cancel) 20 (We thank Barbara Terhal for providing this derivation): Here, the combinatorial factor counts the number of combinations of l errors in k rounds, given an ϵ L chance of error per round. This can be simplified to However, at small k, the decay is dominated by the majority vote, for which ϵ L ∝ (kϵ phys ) (d + 1)/2 . For example, for all the Surface-17 decay curves, we observe a quadratic error rate at small k, as opposed to the linear slope predicted by Eq. (5). In order to correct for this, we shift the above equation in k by a free parameter k 0 , resulting in Eq. (2). This function fits well to data with k ≥ 3 in all plots, and thus allows accurate determination of ϵ L .
The quantumsim simulation package. Quantumsim performs calculations on density matrices utilizing a graphics processing unit in a standard desktop computer. Ancillas are measured at the end of each cycle, and thus not entangled with the rest of the system. As such, it is possible to obtain the effect of the QEC cycle on the system without explicitly representing the density matrix of all 17 qubits simultaneously. The simulation is set up as follows: the density matrix of the nine data qubits is allocated in memory with all qubits initialized to 0 j i. One-qubit and twoqubit gates are applied to the density matrix as completely positive, trace preserving maps represented by Pauli transfer matrices. When a gate involving an ancilla qubit must be performed, the density matrix of the system is dynamically enlarged to include that one ancilla.
Qubit measurements are simulated as projective and following the Born rule, with projection probabilities given by the squared overlap of the input state with the measurement basis states. In order to capture empirical measurement errors, we implement a black-box measurement model (see Error Models) by sandwiching the measurement between idling processes. The measurement projects the system to a product state of the ancilla and the projected sub-block of the density matrix. We can therefore remove the ancilla from the density matrix and only store its state right after projection, and continue the calculation with the partial density matrix of the other qubits. Making use of the specific arrangement of the interactions between ancillas and data qubits in Surface-17, it is possible to apply all operations to the density matrix in such an order (shown in Supplemental Material) that the total size of the density matrix never exceeds 2 10 × 2 10 (nine data qubits plus one ancilla), which allows relatively fast simulation. We emphasize that with the choice of error model in this work, this approach gives the same result as a full simulation on a 17-qubit density matrix. Only the introduction of residual entangling interactions between data and ancilla qubits (which we do not consider in this work) would make the latter necessary. On our hardware (See Supplemental Material), simulating one QEC cycle of Surface-17 with quantumsim takes 25 ms.
We highlight an important advantage of doing density-matrix calculations with quantumsim. We do not perform projective measurements of the data qubits. Instead, after each cycle, we extract the diagonal of the data-qubit density matrix, which represents the probability distribution if a final measurement were performed. We leave the density matrix undisturbed and continue simulation up to k = 20. This is a very useful property of the density-matrix approach, because having a probability distribution of all final readout events greatly reduces sampling noise.
Our measurement model includes a declaration error probability (see Error Models), where the projected state of the ancilla after measurement is not the state reported to the decoder. Before decoding, we thus apply errors to the outcomes of the ancilla projections, and smear the probability distribution of the data qubit measurement. To then determine the fidelity averaged over this probability distribution, we present all 16 possible final Z-type parities to the decoder. This results in 16 different final Pauli updates, allowing us to determine correctness of the decoder for all 512 possible measurement outcomes. These are then averaged over the simulated probability distribution. This produces good results after about 10 4 simulated runs.
A second highlight of quantumsim is the possibility to quantify the suboptimality of the decoder. The fidelity of the logical qubit obtained in these numerical simulations is a combination of the error rates of the physical qubits and the approximations made by the decoder. Full density-matrix simulations make it possible to disentangle these two contributions. Namely, the fidelity is obtained by assigning correctness to each of the 512 possible readouts according to 16 outputs of the decoder, and summing the corresponding probabilities accordingly. If the probabilities are known, it is easy to determine the 16 results that a decoder should output in order to maximize fidelity (i.e., the output of the best-possible decoder). This allows placing a decoder UB F max L on logical fidelity as limited by the physical qubits independent of the decoder. Conversely, it also allows quantifying sub-optimality in the decoder used. In fact, we can make the following reverse statement: if our measurement model did not include a declaration error, then we could use the simulation to find the final density matrix of the system conditioned on a syndrome measurement. From this, the simulation could output exactly the 16 results that give F max L , so that quantumsim could thus be used as a maximum-likelihood decoder. In this situation, F max L would not only be an UB, but indeed the performance of the best-possible decoder. However, as we add the declaration errors after simulation, we can only refer to F max L as the decoder UB.

Error models
We now describe the error model used in the simulations. Our motivation for the development of this error model is to provide a limited number of free parameters to study, whilst remaining as close to known experimental data as possible. As such, we have taken well-established theoretical models as a base, and used experimental tomography to provide fixed parameters for observed noise beyond these models. The parameters of the error model are provided in Supplemental Material.
Density-matrix simulation of small surface TE O'Brien et al.
Idling qubits. While idling for a time τ, a transmon in 1 j i can relax to 0 j i. Furthermore, a transmon in superposition can acquire random quantum phase shifts between 0 j i and 1 j i due to 1/f noise sources (e.g., flux noise) and broadband ones (e.g., photon shot noise 35 and quasiparticle tunneling 36 ). These combined effects can be parametrized by probabilities p 1 = exp(−τ/T 1 ) for relaxation, and p ϕ = exp(−τ/T ϕ ) for pure dephasing. The combined effects of relaxation and pure dephasing lead to decay of the off-diagonal elements of the qubit density matrix. We model dephasing from broadband sources in this way, taking for T ϕ the value extracted from the decay time T 2 of standard echo experiments: We model 1/f sources differently, as discussed below.
Dephasing from photon noise. The dominant broadband dephasing source is the shot noise due to photons in the readout resonator. This dephasing is present whenever the coupled qubit is brought into superposition before the readout resonator has returned to the vacuum state following the last measurement. This leads to an additional, timedependent pure dephasing (rates given in Supplemental Material).
One-qubit Y rotations. We model y-axis rotations as instantaneous rotations sandwiched by idling periods of duration τ g,1Q /2. The errors in the instantaneous gates are modeled from process matrices measured by gate-set tomography 30,37 in a recent experiment. 20 In this experiment, the GST analysis of single-qubit gates also showed that the errors can mostly be attributed to Markovian noise. For simplicity, we thus model these errors as Markovian.
Dephasing of flux-pulsed qubits. During the coherent step, transmons are repeatedly moved in frequency away from their sweetspot using flux pulses, either to implement a C-Z gate or to avoid one. Away from the sweetspot, transmons become first-order sensitive to flux noise, which causes an additional random phase shift. As this noise typically has a 1/f power spectrum, the largest contribution comes from low-frequency components that are essentially static for a single run, but fluctuating between different runs. In our simulation, we approximate the effect of this noise through ensemble averaging, with quasi-static phase error added to a transmon whenever it is flux pulsed. Gaussian phase errors with the variance (calculated in Supplemental Material) are drawn independently for each qubit and for each run.
C-Z gate error. The C-Z gate is achieved by flux pulsing a transmon into the 11 j i $ 02 j i avoided crossing with another, where the 2 denotes the second-excited state of the fluxed transmon. Holding the transmons here for τ g,2Q causes the probability amplitudes of 01 j i and 11 j i to acquire phases. 38 Careful tuning allows the phase ϕ 01 acquired by 01 j i (the single-qubit phase ϕ 1Q ) to be an even multiple of 2π, and the phase ϕ 11 acquired by 11 j i to be π extra. This extra phase acquired by 11 j i is the two-qubit phase ϕ 2Q . Single-qubit and two-qubit phases are affected by flux noise because the qubit is first-order sensitive during the gate. Previously, we discussed the single-qubit phase error. In Supplemental Material, we calculate the corresponding two-qubit phase error δϕ 2Q . Our full (but simplistic) model of the C-Z gate consists of an instantaneous C-Z gate with single-qubit phase error δϕ 1Q and twoqubit phase error δϕ 2Q = δϕ 1Q /2, sandwiched by idling intervals of duration τ g,2Q /2. Measurement. We model qubit measurement with a black-box description using parameters obtained from experiment. This description consists of the eight probabilities for transitions from an input state i j i 2 0 j i; 1 j i f g into pairs m; o j i ð Þof measurement outcome m 2 þ1; À1 f gand final state o 2 0 j i; 1 j i f g . By final state we mean the qubit state following the photondepletion period. Input superposition states in the computational bases are first projected to 0 j i and 1 j i following the Born rule. The probability tree (the butterfly) is then used to obtain an output pair m; o j i ð Þ . These experimental parameters can be described by a six-parameter model (described in detail in Supplemental Material), consisting of periods of enhanced noise before and after a point at which the qubit is perfectly projected, and two probabilities ϵ i j i RO for wrongly declaring the result of this projective measurement. In Supplemental Material, a scheme for measuring these butterfly parameters and mapping them to the six-parameter model is described. In experiment, we find that the readout errors ϵ i j i RO are almost independent of the qubit state i j i, and so we describe them with a single readout error parameter ϵ RO in this work.

Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding author (L.DiCarlo@tudelft.nl) upon request.

Code availability
The computer code that was use to generate these results is available from the corresponding author (L.DiCarlo@tudelft.nl) upon request.