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

TL;DR

This paper uses density-matrix simulations to evaluate the performance of small surface codes under realistic noise, showing they can achieve error rates below fault-tolerance thresholds and potentially surpass quantum memory and computation break-even points.

Contribution

It provides a detailed simulation framework for surface codes with experimental noise parameters and explores optimization strategies and scalability to larger codes.

Findings

Surface-17 achieves error rates below fault-tolerance thresholds.

Feedback strategies improve quantum memory performance.

Surface-49 could surpass the break-even point with current technology.

Abstract

We present a full 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 surface code, and the potential for Surface-17 to perform beyond the break-even point of quantum memory. At state-of-the-art qubit relaxation times and readout speeds, Surface-49 could surpass the break-even point of computation.

Tables

Table 1. Table 1: Standard simulation parameters: Summary of standard times used in all density-matrix simulations, unless otherwise indicated. The two-qubit gate is a conditional phase gate (C-Z). Other error rates and parameters are given in  Suppmaterial .

Parameter	Symbol	Value	Reference
Qubit relaxation time	$T_1$	$30\mu s$	Bultink16
Qubit dephasing time (white noise)	$T_{\varphi }$	$60\mu s$	Asaad16 ; Bultink16
Single-qubit gate time	${\tau }_{g,1Q}$	$20\mathrm{ns}$	Asaad16 ; Bultink16
Two-qubit gate time	${\tau }_{g,2Q}$	$40\mathrm{ns}$	Riste15
Coherent step time	${\tau }_c$	$200\mathrm{ns}$	Versluis16
Measurement time	${\tau }_m$	$300\mathrm{ns}$	Bultink16
Depletion time	${\tau }_d$	$300\mathrm{ns}$	Bultink16
Fast measurement time	${\tau }_m^{\left(\mathrm{fast}\right)}$	$100\mathrm{ns}$	Walter17
Fast depletion time	${\tau }_d^{\left(\mathrm{fast}\right)}$	$100\mathrm{ns}$	Walter17

Table 2. Table 2: Standard parameters of error models used in quantumsim, unless indicated otherwise.

Parameter	Symbol	Value	Reference
In-axis rotation error	$p_{\mathrm{axis}}$	${10}^{-4}$	Rol16
In-plane rotation error	$p_{\mathrm{plane}}$	$5\times {10}^{-4}$	Rol16
$1/f$ flux noise.	$A$	${\left(1\mu {\Phi }_0\right)}^2$	Yan15; Quintana17
Readout infidelity	${ϵ}_{\mathrm{RO}}$	$5\times {10}^{-3}$	Bultink16
Photon relaxation time	$1/\kappa$	$250\mathrm{ns}$	Bultink16
Dispersive shift	$\chi /\pi$	$-2.6\mathrm{MHz}$	Bultink16
photon $\#$ post-measurement	$n_0$	$0.8$ photons	Bultink16

Table 3. Table 3: Measurement butterfly matching a recent characteristic experiment  Bultink16 using a Josephson parametric amplifier  Castellanos-Beltran08 in phase-preserving mode as the front end of the readout amplification chain.

Probability	Value	Probability	Value
${ϵ}_0^{+1,0}$	0.9985	${ϵ}_1^{+1,0}$	0.0050
${ϵ}_0^{+1,1}$	0.0000	${ϵ}_1^{+1,1}$	0.0015
${ϵ}_0^{-1,0}$	0.0015	${ϵ}_1^{-1,0}$	0.0149
${ϵ}_0^{-1,1}$	0.000	${ϵ}_1^{-1,1}$	0.9786

Equations

$${\mathcal{F}}_{\mathrm{phys}}(t)=\frac{1}{6}\left(1+e^{-t/T_1}\right)+\frac{1}{3}\left(1+e^{-t(1/2T_1+1/T_{\varphi })}\right).$$

$${\mathcal{F}}_{\mathrm{L}}[k]=\frac{1}{2}[1+(1-2{ϵ}_{\mathrm{L}}{)}^{k-k_0}].$$

$${ϵ}_{\mathrm{phys}}=-{\tau }_{\mathrm{cycle}}{\frac{d{\mathcal{F}}_{\mathrm{phys}}(t)}{dt}|}_{t=0}=\frac{{\tau }_{\mathrm{cycle}}}{3T_1}+\frac{{\tau }_{\mathrm{cycle}}}{3T_{\varphi }}.$$

$${\mathcal{F}}_{\mathrm{L}}[k]=1-\sum _{l\mathrm{odd}}\left(\genfrac{}{}{0ex}{}{k}{l}\right){ϵ}_{\mathrm{L}}^l(1-{ϵ}_{\mathrm{L}}{)}^{k-l}.$$

$${\mathcal{F}}_{\mathrm{L}}[k]$$

$$=1-\frac{1}{2}\left[(1-{ϵ}_L+{ϵ}_L{)}^k-(1-{ϵ}_L-{ϵ}_L{)}^k\right]$$

$$=\frac{1}{2}[1+(1-2{ϵ}_{\mathrm{L}}{)}^k].$$

$$\frac{1}{T_2}=\frac{1}{T_{\varphi }}+\frac{1}{2T_1}.$$

$${\left(R_{\Lambda }\right)}_{ij}=\frac{1}{2}\text{Tr}\left({\sigma }_i\Lambda {\sigma }_j\right),$$

$$R_{{\Lambda }_{T_1}}$$

$$R_{{\Lambda }_{T_{\varphi }}}$$

$$R_{AP(t)}=R_{{\Lambda }_{T_1}}{R}_{{\Lambda }_{T_{\varphi }}}$$

$$\frac{d\rho }{dt}=-i\frac{B}{2}[{\sigma }_z,\rho ]+\frac{{\Gamma }_d}{2}\mathcal{D}[{\sigma }_z]\rho .$$

$${\alpha }_g{\alpha }_e^{\ast }=\alpha (t_m)\exp \left(-\kappa \left(t-t_m\right)\right)\exp \left(2i\chi \left(t-t_g\right)\right),$$

$$\begin{array}[]{lll}p_{\phi,\mathrm{photon}}&=&\exp\left(-\int_{t_{1}}^{t_{2}}\Gamma_{d}(t)dt\right)\\ &=&\exp\left(\vphantom{{\left[\frac{e^{-\kappa t}}{4\chi^{2}+\kappa^{2}}[-\kappa\sin(2\chi t)-2\chi\cos(2\chi t)]\right]}^{t_{2}-t_{g}}_{t_{1}-t_{g}}}2\chi\alpha(0)\exp(\kappa(t_{m}-t_{g}))\right.\\ &&\times\left.{\left[\frac{e^{-\kappa t}}{4\chi^{2}+\kappa^{2}}[-\kappa\sin(2\chi t)-2\chi\cos(2\chi t)]\right]}^{t_{2}-t_{g}}_{t_{1}-t_{g}}\right).\end{array}$$

$$R_{R_y(\pi /2)}=R_{AP({\tau }_{\mathrm{g},1Q}/2)}{R}_{R_y(\pi /2)}^{\prime }{R}_{\text{dep}}{R}_{AP({\tau }_{\mathrm{g},1Q}/2)},$$

$$R_{\text{dep}}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1-p_{\mathrm{plane}}& 0& 0\\ 0& 0& 1-p_{\mathrm{axis}}& 0\\ 0& 0& 0& 1-p_{\mathrm{plane}}\end{array}\right),$$

$$f_{\mathrm{q}}(t)=(f_{\mathrm{q},\max }+E_{\mathrm{C}})\sqrt{|\cos \left(\pi \Phi (t)/{\Phi }_0\right)|}-E_{\mathrm{C}}$$

$$\frac{\partial f_{\mathrm{q}}}{\partial \Phi }=\frac{-\pi }{2{\Phi }_0}(f_{\mathrm{q}}+E_{\mathrm{C}})\tan \left(\frac{\pi \Phi }{{\Phi }_0}\right).$$

$$\delta \varphi =-2\pi {\tau }_{\mathrm{g},2Q}\frac{\partial f_{\mathrm{q}}}{\partial \Phi }\delta \Phi .$$

$$S_{\Phi }(f)\approx A/f,$$

$$\begin{array}[]{lcl}\delta\Phi_{\mathrm{rms}}&=&{\left(\int_{f_{\min}}^{f_{\max}}S_{\Phi}(f)\,df\right)}^{1/2}\\ &=&A{\left(\ln\left(f_{\max}/f_{\min}\right)\right)}^{1/2}\\ &\approx&4~{}\mu\Phi_{0}.\end{array}$$

$$\delta {\varphi }_{\mathrm{rms}}\approx 0.01\mathrm{rad}.$$

$${\varphi }_{2\mathrm{Q}}={\varphi }_{11}-{\varphi }_{01}=-2\pi {\int }_{t_1}^{t_2}\zeta (t)dt,$$

$$\zeta \approx \beta (\Phi -{\Phi }_{\mathrm{c}})-\sqrt{{\beta }^2{\left(\Phi -{\Phi }_{\mathrm{c}}\right)}^2+{\left(2J/2\pi \right)}^2},$$

$$\beta =\frac{1}{2}\frac{\partial f_{\mathrm{q}}}{\partial \Phi }{|}_{\Phi ={\Phi }_{\mathrm{c}}}.$$

$$\frac{\partial \zeta }{\partial \Phi }{|}_{\Phi ={\Phi }_{\mathrm{c}}}=\beta .$$

$$\delta {\varphi }_{2\mathrm{Q}}\approx -2\pi {\tau }_{\mathrm{g},2Q}\frac{\partial \zeta }{\partial \Phi }{|}_{\Phi ={\Phi }_{\mathrm{c}}}\delta \Phi .$$

$$\begin{array}[]{lcl}\mathrm{P}{(m_{\mathrm{a}}=+1)}_{i}&=&\epsilon_{i}^{+1,0}+\epsilon_{i}^{+1,1},\\ \mathrm{P}{(m_{\mathrm{a}}=+1)}_{i}&=&\epsilon_{i}^{-1,0}+\epsilon_{i}^{-1,1},\\ \mathrm{P}{(m_{\mathrm{b}}=m_{\mathrm{a}}=+1)}_{i}&=&\left(\epsilon_{0}^{+1,0}+\epsilon_{0}^{+1,1}\right)\epsilon_{i}^{+1,0}\\ &&+\left(\epsilon_{1}^{+1,0}+\epsilon_{1}^{+1,1}\right)\epsilon_{i}^{+1,1},\\ \mathrm{P}{(m_{\mathrm{b}}=-m_{\mathrm{a}}=+1)}_{i}&=&\left(\epsilon_{0}^{+1,0}+\epsilon_{0}^{+1,1}\right)\epsilon_{i}^{-1,0}\\ &&+\left(\epsilon_{1}^{+1,0}+\epsilon_{1}^{+1,1}\right)\epsilon_{i}^{-1,1},\\ \mathrm{P}{(-m_{\mathrm{b}}=m_{\mathrm{a}}=+1)}_{i}&=&\left(\epsilon_{0}^{-1,0}+\epsilon_{0}^{-1,1}\right)\epsilon_{i}^{+1,0}\\ &&+\left(\epsilon_{1}^{-1,0}+\epsilon_{1}^{-1,1}\right)\epsilon_{i}^{+1,1},\\ \mathrm{P}{(-m_{\mathrm{b}}=-m_{\mathrm{a}}=+1)}_{i}&=&\left(\epsilon_{0}^{-1,0}+\epsilon_{0}^{-1,1}\right)\epsilon_{i}^{-1,0}\\ &&+\left(\epsilon_{1}^{-1,0}+\epsilon_{1}^{-1,1}\right)\epsilon_{i}^{-1,1}.\end{array}$$

$$\begin{gathered} p_{u,v}\approx \sum _{\text{paths }(e_1,e_2,\dots ,e_n) \\ \text{ between }u\text{ and }v}\prod _{j=1}^n{p}_{e_j}, \end{gathered}$$

$$A_S=A_A+A_A^2+A_A^3+\cdots =\frac{1}{1-A_A}-1,$$

$$A_{S}^{t+1}\approx\left(\begin{array}[]{cc}A_{S}^{t}&C_{S}^{t+1}\\ {(C_{S}^{t+1})}^{T}&{(\mathbf{1}-\partial A_{A}^{t+1})}^{-1}\end{array}\right)$$

Full text

I Introduction

Recent experimental demonstrations of small quantum simulations OMalley16 ; Barends15 ; Langford16 and quantum error correction (QEC) Kelly15 ; Riste15 ; Corcoles15 ; Ofek16 position superconducting circuits for targeting quantum supremacy Boixo16 and quantum fault tolerance Martinis15 , 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 Fowler12 ; Landahl11 ; Yoder16 with different error decoders Tomita14 ; Fowler09 ; Heim16 . 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) Blais04 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 well-approximated 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 Fowler12 ; Tomita14 ; Horsman12 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 Versluis16 and the measured performance of current experimental multi-transmon cQED platforms Bultink16 ; Rol16 ; Asaad16 ; Walter17 . For this purpose, we have developed an open-source density-matrix simulation package named quantumsim 111Please visit https://github.com/brianzi/quantumsim. We use quantumsim to extract the logical error rate per QEC cycle, $\epsilon_{\mathrm{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 -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.

II Results

II.1 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 -5 repetition codes Kelly15 , we first focus on the performance of the logical qubit as a quantum memory. Specifically, we quantify the ability to hold a logical $\left\lvert 0\right\rangle$ state, by initializing this state, holding it for $k\in\left\{1,\ldots,20\right\}$ cycles, performing error correction, and determining a final logical state (see Fig. 6 for details). The logical fidelity $\mathcal{F}_{\mathrm{L}}[k]$ is then given by the probability to match the initial state. We observe identical results when using $\left\lvert 1\right\rangle$ or $\left\lvert\pm\right\rangle=\frac{1}{\sqrt{2}}(\left\lvert 0\right\rangle\pm\left\lvert 1\right\rangle)$ in place of $\left\lvert 0\right\rangle$.

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 Suppmaterial . We focus on the QEC cycle proposed in Versluis16 , which pipelines the execution of $X$- and $Z$-type stabilizer measurements. Each stabilizer measurement consists of three parts: a coherent step (duration $\tau_{\mathrm{c}}=2\tau_{\mathrm{g,1Q}}+4\tau_{\mathrm{g,2Q}}$), measurement ($\tau_{\mathrm{m}}$), and photon depletion from readout resonators ($\tau_{\mathrm{d}}$), making the QEC cycle time $\tau_{\mathrm{cycle}}=\tau_{\mathrm{c}}+\tau_{\mathrm{m}}+\tau_{\mathrm{d}}$.

Simulating this concrete quantum circuit with the listed parameters using quantumsim, we predict $\mathcal{F}_{\mathrm{L}}[k]$ of Surface-$17$ (Fig. 1). We show $\mathcal{F}_{\mathrm{L}}[k]$ for both a homemade MWPM decoder (green, described in Suppmaterial ), and an implementation of the LT decoder of Tomita14 (blue, described in Suppmaterial ). To isolate decoder performance, we can compare the achieved fidelity to an upper bound extractable from the density-matrix simulation (red, described in Sec. IV.1.3). 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:

[TABLE]

The observation of $\mathcal{F}_{\mathrm{L}}[k]>\mathcal{F}_{\mathrm{phys}}(k\tau_{\mathrm{cycle}})$ for large $k$ would constitute a demonstration of QEC beyond the quantum memory break-even point Ofek16 . Equivalently, one can extract a logical error rate $\epsilon_{\mathrm{L}}$ from a best fit to $\mathcal{F}_{\mathrm{L}}[k]$ (as derived in Sec. IV.1.2 as the probability of an odd number of errors occurring),

[TABLE]

Here, $k_{0}$ and $\epsilon_{\mathrm{L}}$ are the parameters to be fit. We compare $\epsilon_{\mathrm{L}}$ to the physical error rate

[TABLE]

We observe $\epsilon_{\mathrm{L}}=1.44\,\%_{\mathrm{c}}$ for the LT decoder, $\epsilon_{\mathrm{L}}=1.07\,\%_{\mathrm{c}}$ for the MWPM decoder, and $\epsilon_{\mathrm{L}}=0.68\,\%_{\mathrm{c}}$ at the decoder upper bound ($\%_{\mathrm{c}}$ = $\%$ per cycle). The latter two fall below $\epsilon_{\mathrm{phys}}=1.33\,\%_{\mathrm{c}}$. Defining the decoder efficiency $\eta_{\mathrm{d}}=\epsilon_{\mathrm{L}}^{\mathrm{(UB)}}/\epsilon_{\mathrm{L}}$, we find $\eta_{\mathrm{d}}^{\mathrm{(LT)}}=0.47$ and $\eta_{\mathrm{d}}^{\mathrm{(MWPM)}}=0.64$.

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$ 222A 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 $\mathcal{F}_{\mathrm{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 look-up table decoder).

With the baseline for current performance established, we next investigate $\epsilon_{\mathrm{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_{\phi}$.

II.2 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 $\epsilon_{\mathrm{L}}$ may be to decrease $\tau_{\mathrm{cycle}}$. Currently, the cycle is dominated by $\tau_{\mathrm{m}}+\tau_{\mathrm{d}}$. At fixed readout power, reducing $\tau_{\mathrm{m}}$ and $\tau_{\mathrm{d}}$ will reduce $\tau_{\mathrm{cycle}}$ at the cost of increased readout infidelity $\epsilon_{\mathrm{RO}}$ (described in Sec. IV.2.6). We explore this trade-off in Fig. 2, using a linear-dispersive readout model FriskKockum12 , keeping $\tau_{\mathrm{m}}=\tau_{\mathrm{d}}$ and assuming no leftover photons. Because of the latter, $\epsilon_{\mathrm{L}}^{\mathrm{(MWPM)}}$ reduces from $1.07\,\%_{\mathrm{c}}$ (Fig. 1) to $0.62\,\%_{\mathrm{c}}$ at $\tau_{\mathrm{m}}=300\,\text{ns}$. The minimum $\epsilon_{\mathrm{L}}^{\mathrm{(MWPM)}}=0.55\,\%_{\mathrm{c}}$ is achieved at around $\tau_{\mathrm{m}}=260~{}\mathrm{ns}$. This is perhaps counterintuitive, as $\epsilon_{\mathrm{phys}}$ reduces only $0.13\,\%_{\mathrm{c}}$ while $\epsilon_{\mathrm{RO}}$ increases $0.5\,\%$. However, it reflects the different sensitivity of the code to different types of errors. Indeed, $\epsilon_{\mathrm{L}}^{\mathrm{(MWPM)}}$ is smaller for $\tau_{\mathrm{m}}=200~{}\mathrm{ns}$ than for $\tau_{\mathrm{m}}=300~{}\mathrm{ns}$, even though $\epsilon_{\mathrm{RO}}$ increases to $5\,\%$. It is interesting to note that the optimal $\tau_{\mathrm{m}}$ for quantum memory, which minimizes logical error per unit time, rather than per cycle, is $\tau_{\mathrm{m}}=280\,\text{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 $\epsilon_{\mathrm{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}(\pi/2)$ gate at the end of the coherent step with a $R_{y}(-\pi/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 upper bound. We observe $\epsilon_{\mathrm{L}}$ improve only $0.05\,\%_{\mathrm{c}}$ in both cases. Future experiments might opt not to pursue these small gains in view of the technical challenges added by feedback control.

II.3 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 $\tau_{\mathrm{cycle}}=800~{}\mathrm{ns}$, and adjust $T_{1}$ values with both $T_{\phi}=2T_{1}$ (common in experiment) and $T_{\phi}=\infty$ (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 $\epsilon_{\mathrm{L}}$ and $\epsilon_{\mathrm{phys}}$ over the $T_{1}$ range covered. For the MWPM decoder (upper bound) and $T_{\phi}=2T_{1}$, the memory figure of merit $\gamma_{\mathrm{m}}=\epsilon_{\mathrm{phys}}/\epsilon_{\mathrm{L}}$ increases from $1.3$ $(2)$ at $T_{1}=30~{}\mu\mathrm{s}$ to $2$ $(5)$ at $100~{}\mu\mathrm{s}$. Completely eliminating white-noise dephasing will increase $\gamma_{\mathrm{m}}$ by $10\%$ with MWPM and $30\%$ at the upper bound.

A key question for any QEC code is how $\epsilon_{\mathrm{L}}$ scales with code distance $d$. Computing power limitations preclude similar density-matrix 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 Suppmaterial for details). Figure 5 shows the lowest-order approximation to the logical error rates of Surface-$17$ and -$49$ over a range of $T_{1}=T_{\phi}/2$. Comparing the Surface-$17$ lowest-order approximation to the quantumsim result shows good agreement and validates the approximation. We observe a lower $\epsilon_{\mathrm{L}}$ for Surface-49 than for -17, indicating quantum fault tolerance over the $T_{1}$ range covered. The fault-tolerance figure of merit defined in Martinis15 , $\Lambda_{\mathrm{t}}=\epsilon_{\mathrm{L}}^{(17)}/\epsilon_{\mathrm{L}}^{(49)}$, increases from $2$ to $4$ as $T_{1}$ grows from $30$ to $100~{}\mu\mathrm{s}$.

As a rough metric of computational performance, we offer to compare $\epsilon_{\mathrm{L}}$ (per cycle) to the error accrued by a physical qubit idling over $\tau_{\mathrm{g,1Q}}$. We define a metric for computation performance, $\gamma_{\mathrm{c}}=(\epsilon_{\mathrm{phys}}\tau_{\mathrm{g,1Q}})/(\epsilon_{\mathrm{L}}\tau_{\mathrm{cycle}})$ and $\gamma_{\mathrm{c}}=1$ as a computational break-even point. Clearly, using the QEC cycle parameters of Table 1 and even with $T_{1}$ improvements, neither Surface-17 nor -49 can break-even computationally. However, including the readout acceleration recently demonstrated in Walter17 , which allows $\tau_{\mathrm{m}}=\tau_{\mathrm{d}}=100~{}\mathrm{ns}$ and $\tau_{\mathrm{cycle}}=400~{}\mathrm{ns}$, Surface-49 can cross $\gamma_{\mathrm{c}}=1$ by $T_{1}=40~{}\mu\mathrm{s}$. In view of first reports of $T_{1}$ up to $80~{}\mu\mathrm{s}$ emerging for planar transmons Paik16 ; Gustavsson16 , this important milestone may be within grasp.

III Discussion

III.1 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 Horsman12 ; Yoder16 ), 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.

III.2 Decoder performance

A practical question facing quantum error correction is how best to balance the trade-off between decoder complexity and performance. Past proposals for surface-code computation via lattice surgery Horsman12 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 Terhal15 . 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 $\tau_{\mathrm{cycle}}<1\,\mu\mathrm{s}$. This short time makes the use of computationally intensive decoding schemes difficult, even if they provide lower $\epsilon_{\mathrm{L}}$.

The leading strategy for decoding the surface code is MWPM using the blossom algorithm of Edmonds Fowler12 ; Fowler09 ; Fowler14 . Although this algorithm is challenging to implement, it scales linearly in code distance Fowler14 . 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 Suppmaterial ) is whether these weights can be calculated on the fly, or must be precalculated and stored. On-the-fly 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, Suppmaterial shows that MWPM is nearly perfect in the absence of $Y$ errors. The decoder efficiency $\eta_{\mathrm{d}}$ may significantly increase by extending MWPM to account for correlations between detected $X$ and $Z$ errors originating from $Y$ errors Delfosse14 ; Fowler13b .

If computational limitations preclude a MWPM decoder from keeping up with $\tau_{\mathrm{cycle}}$, the look-up table decoder may provide a straightforward solution for Surface-17. However, at current physical performance, the $\eta_{\mathrm{d}}$ reduction will make Surface-17 barely miss memory break-even (Fig. 1). Furthermore, memory requirements make look-up table 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.

III.3 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 the $T_{1}$ noise, but a trace distance $100$ times smaller. As the flux noise increases $\epsilon_{\mathrm{L}}$ by only $0.01\,\%_{\mathrm{c}}$, it appears $\epsilon_{\mathrm{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 $\epsilon_{\mathrm{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 ($\epsilon_{\mathrm{RO}}$) or the ancilla error rate ($\epsilon_{\mathrm{anc}}$) is below the data qubit ($\epsilon_{\mathrm{phys}}$) error rate, $\epsilon_{\mathrm{L}}\propto(\epsilon_{\mathrm{anc}}+\epsilon_{\mathrm{RO}})\epsilon_{\mathrm{phys}}^{(d-1)/2}$. However, if $\epsilon_{\mathrm{RO}}$ ($\epsilon_{\mathrm{anc}}$) $>\epsilon_{\mathrm{phys}}$, $\epsilon_{\mathrm{L}}$ becomes independent of $\epsilon_{\mathrm{RO}}$ ($\epsilon_{\mathrm{anc}}$), to lowest order. This can be seen in Fig. 2, where the error rate is almost constant as $\epsilon_{\mathrm{RO}}$ exponentially increases. This approximation breaks down with large enough $\epsilon_{\mathrm{anc}}$ and $\epsilon_{\mathrm{RO}}$, but presents a counterintuitive point for experimental design; $\epsilon_{\mathrm{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 $\tau_{\mathrm{cycle}}=720~{}\mathrm{ns}$ to minimize $\epsilon_{\mathrm{L}}$ would lead to a $1.09\%$ error rate for the two qubits combined. However, shortening the $\tau_{\mathrm{cycle}}$ of $q_{c}$ reduces the time over which $q_{m}$ decays. If $q_{c}$ operates at $\tau_{\mathrm{cycle}}=600~{}\mathrm{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.

III.4 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 Blume13 ). 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 to $\sim 0.3\%$ Barends14 , and per single-qubit gate to $\sim 0.001\%$ Chen16 . Schemes have also been developed to reduce the accumulation of leakage Fowler13 . 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}\approx 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.

IV Methods

IV.1 Simulated experimental procedure

IV.1.1 Surface-$17$ basics

A QEC code can be defined by listing the data qubits and the stabilizer measurements that are repeatedly performed upon them GottesmanPhD . In this way, Surface-$17$ is defined by a $3\times 3$ grid of data qubits $\{D_{0},\ldots 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$- and $Z$-type parity operators $X_{2}X_{1}$, $Z_{3}Z_{0}$, $X_{4}X_{3}X_{1}X_{0}$, $Z_{5}Z_{4}Z_{2}Z_{1}$, $Z_{7}Z_{6}Z_{4}Z_{3}$, $X_{8}X_{7}X_{5}X_{4}$, $Z_{8}Z_{5}$, 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 lattices, followed by ancilla measurements [Fig. 6(a)]. 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 Versluis16 , for which we give a schematic description in Fig. 6(b) (see Suppmaterial 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 (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 Fowler12 . The rate at which these logical errors arise is the main focus of this paper.

IV.1.2 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 $\left\lvert 0\right\rangle$ and perform $k=1,2,\ldots,20$ QEC cycles [Fig. 6(b)]. 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 $\left\lvert 0\right\rangle$ state Fowler12 . 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. 6(c)] are processed using a classical decoding algorithm. The decoder computes a Pauli update after each QEC cycle [Fig. 6(d)]. 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 $\left\lvert 0\right\rangle$ 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. 6(e)]. 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 $\mathcal{F}_{\mathrm{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 $\epsilon_{\mathrm{L}}$. The fidelity $\mathcal{F}_{\mathrm{L}}[k]$ is obtained by counting the probability of an odd number of errors having occurred in total (as two $\sigma_{x}$ errors cancel) Rol16 ; Terhalcomment :

[TABLE]

Here, the combinatorial factor counts the number of combinations of $l$ errors in $k$ rounds, given an $\epsilon_{\mathrm{L}}$ chance of error per round. This can be simplified to

[TABLE]

However, at small $k$, the decay is dominated by the majority vote, for which $\epsilon_{\mathrm{L}}\propto(k\epsilon_{\mathrm{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\geq 3$ in all plots, and thus allows accurate determination of $\epsilon_{\mathrm{L}}$.

IV.1.3 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 $\left\lvert 0\right\rangle$. One- and two-qubit 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 (Sec. IV.2.6) 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 Suppmaterial ) that the total size of the density matrix never exceeds $2^{10}\times 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 Suppmaterial ), simulating one QEC cycle of Surface-17 with quantumsim takes $25~{}\mathrm{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 Sec. IV.2.6), 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 $\sim 10^{4}$ simulated runs.

A second highlight of quantumsim is the possibility to quantify the sub-optimality 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 upper bound $\mathcal{F}_{\mathrm{L}}^{\max}$ 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 $\mathcal{F}_{\mathrm{L}}^{\max}$, so that quantumsim could thus be used as a maximum-likelihood decoder. In this situation, $\mathcal{F}_{\mathrm{L}}^{\max}$ would not only be an upper bound, but indeed the performance of the best-possible decoder. However, as we add the declaration errors after simulation, we can only refer to $\mathcal{F}_{\mathrm{L}}^{\max}$ as the decoder upper bound.

IV.2 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 Suppmaterial .

IV.2.1 Idling qubits

While idling for a time $\tau$, a transmon in $\left\lvert 1\right\rangle$ can relax to $\left\lvert 0\right\rangle$. Furthermore, a transmon in superposition can acquire random quantum phase shifts between $\left\lvert 0\right\rangle$ and $\left\lvert 1\right\rangle$ due to $1/f$ noise sources (e.g., flux noise) and broadband ones (e.g., photon shot noise Sears12 and quasiparticle tunneling Riste13 ). These combined effects can be parametrized by probabilities $p_{1}=\exp(-\tau/T_{1})$ for relaxation, and $p_{\phi}=\exp(-\tau/T_{\phi})$ 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_{\phi}$ the value extracted from the decay time $T_{2}$ of standard echo experiments:

[TABLE]

We model $1/f$ sources differently, as discussed below.

IV.2.2 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, time-dependent pure dephasing (rates given in Suppmaterial ).

IV.2.3 One-qubit Y rotations

We model $y$-axis rotations as instantaneous rotations sandwiched by idling periods of duration $\tau_{\mathrm{g,1Q}}/2$. The errors in the instantaneous gates are modeled from process matrices measured by gate-set tomography Blume13 ; Blume16 in a recent experiment Rol16 . 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.

IV.2.4 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 Suppmaterial ) are drawn independently for each qubit and for each run.

IV.2.5 C-Z gate error

The C-Z gate is achieved by flux pulsing a transmon into the $\left\lvert 11\right\rangle\leftrightarrow\left\lvert 02\right\rangle$ avoided crossing with another, where the $2$ denotes the second-excited state of the fluxed transmon. Holding the transmons here for $\tau_{\mathrm{g,2Q}}$ causes the probability amplitudes of $\left\lvert 01\right\rangle$ and $\left\lvert 11\right\rangle$ to acquire phases DiCarlo09 . Careful tuning allows the phase $\phi_{01}$ acquired by $\left\lvert 01\right\rangle$ (the single-qubit phase $\phi_{1\mathrm{Q}}$) to be an even multiple of $2\pi$, and the phase $\phi_{11}$ acquired by $\left\lvert 11\right\rangle$ to be $\pi$ extra. This extra phase acquired by $\left\lvert 11\right\rangle$ is the two-qubit phase $\phi_{2\mathrm{Q}}$. Single- 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 Suppmaterial , we calculate the corresponding two-qubit phase error $\delta\phi_{2\mathrm{Q}}$. Our full (but simplistic) model of the C-Z gate consists of an instantaneous C-Z gate with single-qubit phase error $\delta\phi_{1\mathrm{Q}}$ and two-qubit phase error $\delta\phi_{2\mathrm{Q}}=\delta\phi_{1\mathrm{Q}}/2$, sandwiched by idling intervals of duration $\tau_{\mathrm{g,2Q}}/2$.

IV.2.6 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 $\left\lvert i\right\rangle\in\{\left\lvert 0\right\rangle,\left\lvert 1\right\rangle\}$ into pairs ($m$,$\left\lvert o\right\rangle$) of measurement outcome $m\in\{+1,-1\}$ and final state $\left\lvert o\right\rangle\in\{\left\lvert 0\right\rangle,\left\lvert 1\right\rangle\}$. By final state we mean the qubit state following the photon-depletion period. Input superposition states in the computational bases are first projected to $\left\lvert 0\right\rangle$ and $\left\lvert 1\right\rangle$ following the Born rule. The probability tree (the butterfly) is then used to obtain an output pair $\left(m,\left\lvert o\right\rangle\right)$. These experimental parameters can be described by a six-parameter model (described in detail in Suppmaterial ), consisting of periods of enhanced noise before and after a point at which the qubit is perfectly projected, and two probabilities $\epsilon_{\mathrm{RO}}^{\left\lvert i\right\rangle}$ for wrongly declaring the result of this projective measurement. In Suppmaterial , 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 $\epsilon_{\mathrm{RO}}^{\left\lvert i\right\rangle}$ are almost independent of the qubit state $\left\lvert i\right\rangle$, and so we describe them with a single readout error parameter $\epsilon_{\mathrm{RO}}$ in this work.

Acknowledgements.

We thank C. C. Bultink, M. A. Rol, B. Criger, X. Fu, S. Poletto, R. Versluis, P. Baireuther, D. DiVincenzo, B. Terhal, and C.W.J. Beenakker for useful discussions. This research is supported by the Foundation for Fundamental Research on Matter (FOM), the Netherlands Organization for Scientific Research (NWO/OCW), an ERC Synergy Grant, and by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office grant W911NF-16-1-0071. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon.

Appendix A Full circuit diagram for Surface-17 implementation

The quantum circuit Versluis16 (Fig. 7) consists of $R_{y}(\pi/2)$ (“$+$”) and $R_{y}(-\pi/2)$ (“$-$”) rotations, C-Z gates, and ancilla measurements. The coherent steps of the $X$ and $Z$ ancillas are pipelined (shifted in time with respect to each other) to prevent transmon-transmon avoided crossings. As long as $\tau_{\mathrm{m}}+\tau_{\mathrm{d}}\geq\tau_{\mathrm{c}}$, no time is lost due to this separation.

In a simulation of the given circuit, gates on different qubits commute and may be applied to the density matrix in any order, regardless of the times at which they are performed in an experiment. As described in Sec. IV A 3, by simulating gates in a specific order (Fig. 8), one can ensure that only one ancilla is ancilla is entangled with the data qubits at any point in the simulation. This allows a reduction in the maximum size of the density matrix from $2^{17}\times 2^{17}$ to $2^{10}\times 2^{10}$.

Appendix B Parameters of error models

This appendix provides mathematical details of the sources of error described in the main text. Standard values for the parameters used throughout the text are given in Table 2.

In the quantumsim module, all gates are applied in the Pauli transfer matrix representation Chow12. These are given in the form

[TABLE]

where matrices $\sigma_{i}$ are the Pauli operators: $\sigma_{0}=I$, $\sigma_{1}=X$, $\sigma_{2}=Y$ and $\sigma_{3}=Z$.

B.1 Qubit idling

Idling qubits are described by the amplitude-phase damping model Nielsen00, corresponding to the transfer matrices

[TABLE]

Idling for a duration $t$ is thus described by

[TABLE]

with $p_{1}=1-e^{-t/T_{1}}$ and $p_{\phi}=1-e^{-t/T_{\phi}}$.

B.2 Photon decay

In the presence of photons in a readout resonator, the coupled qubit is affected according to the effective stochastic master equation FriskKockum12 :

[TABLE]

Here, $\rho$ is the qubit density matrix, $\mathcal{D}[X]$ is the Lindblad operator $\mathcal{D}[X]\rho=X\rho X^{{\dagger}}-\frac{1}{2}X^{\dagger}X\rho-\frac{1}{2}\rho X^{\dagger}X$, $B=2\chi\text{Re}(\alpha_{g}\alpha_{e}^{*})$ is the measurement-induced detuning (Stark shift), and $\Gamma_{d}=2\chi\text{Im}(\alpha_{g}\alpha_{e}^{*})$ is the measurement-induced dephasing, with $\alpha_{i}$ the qubit-state-dependent photon field in the resonator and $2\chi$ the qubit frequency shift per photon. At time $t-t_{g}$ after the qubit superposition is created,

[TABLE]

with $t-t_{m}$ the time since the end of measurement excitation pulse. Integrating over the interval $\left[t_{1},t_{2}\right]$ gives a dephasing term with coefficient

[TABLE]

This dephasing is then implemented via the same Pauli transfer matrix as (17).

B.3 Single-qubit $R_{y}(\pi/2)$ rotations

Single-qubit rotations are modeled by sandwiching an instantaneous Pauli transfer matrix, representing the rotation, with periods of duration $\tau_{\mathrm{g,1Q}}/2$ of amplitude and phase damping. This allows to model the gate for different $T_{1}$ and $T_{\phi}$. However, comparison of this model with Pauli transfer matrices obtained from gate-set tomography experiments shows that actual gates are more accurately described when adding a phenomenological depolarizing noise to the instantaneous part. In the Bloch sphere, this decay corresponds to shrinking toward the origin, with factor $1-p_{\mathrm{axis}}$ along the $y$ axis and $1-p_{\mathrm{plane}}$ along the $x$- and $z$-axes. We thus model

[TABLE]

where

[TABLE]

and $R^{\prime}_{R_{y}(\pi/2)}$ is the Pauli transfer matrix describing a perfect $\pi/2$ rotation around the $y$ axis.

B.4 Flux noise

Shifting the transmon from its sweetspot $f_{\mathrm{q,max}}$ to a lower frequency

[TABLE]

makes it first-order sensitive to flux noise, with sensitivity

[TABLE]

Here, $\Phi$ is the flux bias and $\Phi_{0}=h/2e$ is the flux quantum. For a deviation of $\delta\Phi$, the pulsed transmon incurs a phase error

[TABLE]

Flux noise has a characteristic (single-sided) spectral density

[TABLE]

where $A\approx{(1~{}\mu\Phi_{0})}^{2}$ with $f$ in $\mathrm{Hz}$. We model this noise as quasi-static over the duration ($1/f_{\min}\sim 20~{}\mu\mathrm{s}$, or 20 QEC cycles) of individual runs, but fluctuating between subsequent runs ($1/f_{\max}\sim 20~{}\sec$, or $10^{5}$ runs at $200~{}\mu\mathrm{s}$ intervals). The root-mean-square (rms) fluctuations of flux are therefore

[TABLE]

For our quantum circuit based on Versluis16 , we estimate the corresponding rms phase error induced in a pulsed transmon to be

[TABLE]

B.5 C-Z gates

We now focus on the two-qubit phase error. For an adiabatic gate,

[TABLE]

with $t_{1}$ and $t_{2}=t_{1}+\tau_{\mathrm{g,2Q}}$ the start and end of the gate and $\zeta$ the time-dependent frequency deviation of the lower branch of the $\left\lvert 11\right\rangle\leftrightarrow\left\lvert 02\right\rangle$ avoided crossing from the sum of frequencies for $\left\lvert 01\right\rangle$ and $\left\lvert 10\right\rangle$. Near the flux center $\Phi_{\mathrm{c}}$ of the $\left\lvert 11\right\rangle-\left\lvert 02\right\rangle$ avoided crossing,

[TABLE]

where $2J/2\pi\sim 50~{}\mathrm{MHz}$ is the minimum splitting between $\left\lvert 11\right\rangle$ and $\left\lvert 02\right\rangle$, and

[TABLE]

Differentiating with respect to $\Phi$ at $\Phi_{\mathrm{c}}$ gives

[TABLE]

To estimate the $\delta\phi_{2\mathrm{Q}}$ error, we make the following simplification: we replace the exact trajectory created by the flux pulse by a shift to $\Phi=\Phi_{\mathrm{c}}+\delta\Phi$ with duration $\tau_{\mathrm{g,2Q}}$. For a deviation of $\delta\Phi$,

[TABLE]

Note that this two-qubit phase error is correlated with the single-qubit phase error on the fluxed transmon. The former is smaller by a factor $\approx 2$.

B.6 Measurement

The probabilities $\epsilon_{i}^{m,o}$ are calibrated using the statistics of outcomes in back-to-back measurements ($a$ followed by $b$) with the qubit initialized in $\left\lvert i\right\rangle$.

[TABLE]

We obtain the six free parameters of the black-box description from these 12 equations, using experimental values on the left-hand side Riste12. Table 3 shows the values used, achieved in a recent experiment Bultink16 .

For the simulation, we reproduce this behaviour of the measurement process by a model with several steps. The qubit undergoes dephasing, followed by periods of decay or excitation between which the measurement result is sampled. This measurement result is further subject to a state-dependent declaration error $\epsilon_{\mathrm{RO}}$ before reported to the decoder (see Fig.9). The six parameters of this model are in a one-to-one correspondence with the butterfly parameters described above, and can be mapped by solving the corresponding system of equations. The experimental results in Tab.3 are very well explained by assuming unmodified amplitude-phase damping (withe zero excitation probabilities) during the measurement period, and an outcome-independent declaration error of $\epsilon_{\mathrm{RO}}=\epsilon_{\mathrm{RO}}^{1}=\epsilon_{\mathrm{RO}}^{0}=0.15\%$. We use this result to extrapolate measurement performance to different values of $T_{1}$.

Reduction of measurement time is expected to reduce assignment fidelity. For the results presented in Fig. 2, we do not rely on experimental results, but assume a simplified model for measurement, following Ref. FriskKockum12, . A constant drive pulse of amplitude $\epsilon$ and tuned to the bare resonator frequency, $\Delta_{r}=0$, excites the readout resonator for time $\tau_{\mathrm{m}}$. The dynamics of the resonator is dependent on the transmon state (we approximate linear behavior), and the transmitted signal is amplified and detected in a homodyne measurement as a noisy transient. This transient is processed by a linear classifier, which declares the measurement outcome. For resonator depletion, we use a two-step clearing pulse with amplitude $\epsilon_{c1}$ and $\epsilon_{c2}$, each active for $\tau_{\mathrm{d}}/2$ and chosen (by numerical minimization) so that, at the end of the depletion pulse, the transients for both transmon states return to zero. While the resonator dynamics is easily found if the transmon is in the ground state, amplitude damping of the transmon in the excited state leads to non-deterministic behavior. We thus numerically obtain an ensemble of noisy transients for each input qubit state, and optimize the decision boundary of the linear classifier for this ensemble. Generating a second verification ensemble, the “butterfly” of the measurement setup is estimated.

The dynamics of the resonator is determined by the resonator linewidth $\kappa$ as well as the dispersive shift $\chi$. We chose the parameters of the setup used in Bultink16 , $1/\kappa=250\,\text{ns}$ and $\chi/\pi=-2.6\,\text{MHz}$. The signal-to-noise ratio of the detected transient is reduced by the quantum efficiency $\eta=12.5\%$. The driving strength $\epsilon$ is chosen to approximate the “butterfly” used in most of the main text, and corresponds to a steady-state average photon population of about $\bar{n}=15$. We then keep $\epsilon$ constant while changing the measurement time, keeping $\tau_{\mathrm{m}}=\tau_{\mathrm{d}}$, to obtain the butterflies used in the density matrix simulation. We ignore effects leading to measurement-induced mixing and non-linearity of the readout resonator. Finally, since these simulations do not allow to make a realistic prediction about residual photon numbers achievable in experiments, we ignore this effect when using these results.

Appendix C Effect of over-rotations and two-qubit phase noise on logical error rate

In this section we provide additional numerical data showing the effect of some common noise sources on the logical error rate. In Fig. 10 we show the effect of a coherent over-rotation, whereby the $R^{\prime}_{Y}(\pi/2)$ operator in Eq. 19 is replaced by $R^{\prime}_{Y}(\pi/2+\delta\phi)$. This can be caused by inaccurate calibration of the flux pulse used to perform the gate. In Fig. 11 we show the effect of an increase in the two-qubit flux noise $\delta\phi_{\mathrm{rms}}$ as described in Sec. B.4.

Appendix D Calculation of decoder upper bound

We provide a detailed description how the decoder upper bound is obtained from the simulation results. As described in the main text, after each cycle of simulation, the diagonal of the reduced density matrix of the data qubits in the Z basis is stored. It contains the probability distribution for the $2^{9}=512$ different possible measurement outcomes of the data qubits. In the quantum memory experiment described in the main text, each of these outcomes are passed to the decoder, which then declares a logical measurement outcome.

It is evident that any decoder must declare opposite logical outcomes if two of the 512 possible measurements m and m’ are related by the application of a logical X operator. Thus, any decoder can give the correct result only for half of the measurement outcomes. Subject to this constraint, we can find the set of $256$ declarations which maximize the probability that the declaration is correct. It immediately follows that no decoder can achieve a declaration fidelity larger than this maximal probability. We thus refer to it as the decoder upper bound.

In practice, the upper bound is found according to the following approach. Since declarations are opposite if two outcomes differ by a logical X operator, they must be equal if they differ by the application of one or more X stabilizers (applying two different logical X operators amounts to the application of a product of X stabilizers). We thus group the outcomes in 32 cosets which are related by the application of X-stabilizers. (There are 4 X-stabilizers in Surface-17, so there are $512/2^{4}=32$ cosets). For outcomes from the same coset, the declaration from a decoder must be the same. We obtain the probability of a final measurement falling within each coset by summing the probabilities from the density matrix diagonal. We further group the 32 cosets to 16 pairs, which differ by the application of a logical operator. The upper bound is then obtained by selecting the more probable coset from each pair and summing the corresponding probabilities. This upper bound can also be interpreted as the internal decoherence of the logical qubit: it represents the maximal overlap of the final state with the initial state, under any possible correction of errors.

We finally emphasize that the this upper bound can be found only because we have access to the complete probability distribution of outcomes (for a given result of syndrome measurements), a major advantage of the density matrix simulation. However, we do not expect that any decoder can actually achieve this upper bound: This is because we add syndrome measurement events independently after the situation, which will decrease the logical error rate further.

Appendix E Hardware requirements of simulation

The simulations are performed using the quantumsim package quantumsim_website , which were developed by the authors for this work. The package is accelerated by performing the density matrix manipulations on a GPU (graphics card). The simulations for this work were performed on a NVidia Tesla K40 GPU, on which we observed runtimes of about 0.5 seconds for the simulation of a run of k=20 cycles (25 ms per QEC cycle). We also had the opportunity to test the software on a more modern GPU (NVidia Tesla P100), observing about 15 ms per cycle, and on a consumer-grade GPU (NVidia Quadro M2000), observing about 40 ms per cycle. By comparison, the CPU is mostly idle during the simulation, except for handling of input and output. The memory requirements are modest for both CPU and GPU RAM. They are dominated by the storage of the density matrices and amount to a few ten megabytes.

Appendix F Homemade MWPM decoder with asymmetric weight calculation

Every QEC code requires a decoder to track the most likely errors consistent with a given set of stabilizer measurements. The MWPM decoder has gained popularity since it was shown to have threshold values above $1\%$ Fowler09 . The motivation behind MWPM is that single $X$ or $Z$ errors on data qubits in the bulk of a surface-code fabric cause changes of two stabilizers in the code. These signals can then be considered vertices on a graph, with the error the edge connecting them. Errors in measurement, or errors on a single ancilla qubit, behave as changes in the stabilizer that are separated in time. Multiple errors that would join the same vertices create longer paths in the graph, of which an experiment only records the endpoints. Thus, the problem becomes that of finding the most likely set of generating errors given the error signals that mark their ends. This is made slightly simpler, as in the surface code any chain of errors that forms a closed loop does not change the logical state. This implies that all paths that connect two points are equivalent, and can be considered together. The problem then is to join error signals, either in pairs, or to a ‘boundary’ vertex. The latter corresponds to errors on data qubits at the boundary, which belong to only one $X$ or $Z$ stabilizer. This pairing $P$ should be chosen as the most likely combination of single-qubit errors that could generate the measured error signals. This has then been reduced to the problem of minimum-weight perfect matching on a graph, which can be solved in polynomial time by the blossom algorithm Edmonds65; Fowler12 .

The MWPM decoder we use differs from previous methods by its weight calculation. As part of the decoding process, it is required to calculate to some degree of accuracy Fowler12b the probability $p_{e_{1},e_{2}}$ of two measured error signals $e_{1}$ and $e_{2}$ being connected by a chain of individual logical errors. This is then converted to a weight $w_{e_{1},e_{2}}=-\log(p_{e_{1},e_{2}})$, which form the input to the blossom algorithm of Edmonds to find the most likely matching of error signals Edmonds65; Fowler12 . An exact calculation of $p_{e_{1},e_{2}}$ requires a sum over all such chains between $e_{1}$ and $e_{2}$ that do not cross the boundary (these are equivalent modulo stabilizer operators that do not change the logical state). In this appendix we detail a method of computing this sum, and approximations to make it viable within the runtime of the experiment.

Let us define the ancilla graph $\mathcal{G}_{A}=(V_{A},E_{A})$ containing a vertex $v\in V_{A}$ for every ancilla measurement, and an edge $e\in E_{A}$ connecting $v,u\in V_{A}$ if a single component (gate, single-qubit rest period, or faulty measurement) in the simulation can cause the $u$ and $v$ measurements to return an error. We include a special ‘boundary’ vertex $v_{B}$, to which we connect another vertex $v$ if single components can cause errors on $v$ alone. Then, to each edge $e$ we associate a probability $p_{e}$, being the sum of the probabilities of each component causing this error signal. These error rates can be obtained directly from quantumsim, by cutting the circuit at each C-Z gate and measuring the decay of single qubits between. Then, for a given experiment with given syndrome measurements, let us define the syndrome graph $\mathcal{G}_{S}=(V_{S},E_{S})$ containing a vertex $v\in V_{S}$ for each syndrome measurement that records an error, and an edge $\lambda_{u,v}\in E_{S}$ connecting $u,v\in V_{S}$ if $u$ and $v$ are either both $X$ ancilla qubits or both $Z$ ancilla qubits. To each edge $\lambda_{u,v}$ we associate a probability $p_{u,v}$ given by the sum of the probabilities of a chain of errors causing error signals solely on $u$ and $v$.

If we assume that single-qubit errors are uncorrelated, we have to lowest order

[TABLE]

Let $A_{A}$ be the adjacency matrix on $\mathcal{G}_{A}$ weighted by the probabilities $p_{e}$ (i.e., $(A_{A})_{u,v}=p_{e}$ with $e$ connecting $u$ and $v$), and $A_{S}$ the same for $\mathcal{G}_{S}$. Then, the above becomes

[TABLE]

noting that $A_{S}$ contains a subset of the indices that are used to construct $A_{A}$.

The boundary must be treated specially in the above calculation. For the purposes of the surface code, the boundary can be described as a single vertex which has no limit on the number of other vertices it may pair to Fowler12 . For the purposes of weight calculation, any path that passes through the boundary is already counted by pairing both end vertices to the boundary. This can be treated by making $\mathcal{G}_{A}$ directed, and breaking the symmetry $A_{A}^{T}=A_{A}$. In particular, either ${(A_{A})}_{v_{B},u}=0$ for all $u$ or ${(A_{A})}_{u,v_{B}}=0$ for all $u$.

The above calculation requires inversion of a $N_{\mathrm{mat}}\times N_{\mathrm{mat}}$ matrix, with $N_{\mathrm{mat}}$ the total number of ancilla measurements per experiment. Furthermore, as ancilla error rates depend upon the previous ancilla state, elements in $A_{A}$ are not completely known until the previous cycle. This implies that in an actual computation with runtime decoding, this inversion would need to be completed within a few microseconds (with a transmon-cQED architecture), which is practically unfeasible. We suggest two approximations that can be made to shorten the decoding time. The first is to average all errors over the ancilla population, ignoring any asymmetry in the system. The adjacency matrix is now the same for any experiment, and can be precalculated and stored as a look-up table for the run-time decoder. We call this the decoder with symmetrized weights. The size of such a look-up table scales poorly with the number of qubits and the number of cycles. However, $(A_{S})_{u,v}$ is approximately invariant under simultaneous translation of $u$ and $v$ (excluding boundary effects). This implies that a precalculated $A_{S}$ can be vastly compressed, making this method feasible.

The second approximation to the full $A_{S}$ calculation is to perform it iteratively. We divide our graph $\mathcal{G}_{A}$ ($\mathcal{G}_{S}$) by time steps; let $\mathcal{G}_{A}^{t}$ ($\mathcal{G}_{S}^{t}$) be the subgraph of $\mathcal{G}_{A}$ ($\mathcal{G}_{S}$) containing only ancillas measured before time step $t$, and let $\partial\mathcal{G}_{A}^{t}$ ($\partial\mathcal{G}_{S}^{t}$) be the subgraph of $\mathcal{G}_{A}$ ($\mathcal{G}_{S}$) containing only ancillas measured during time step $t$. Then, if we assume we have an approximation to the matrix $A_{S}^{t}$ (being the adjacency matrix of $\mathcal{G}_{S}^{t}$), we can approximate

[TABLE]

to lowest order in physical errors. Here, $\partial A_{A}^{t+1}$ is the weighted adjacency matrix on $\partial\mathcal{G}_{A}^{t+1}$, and the coupling matrix $C_{S}^{t+1}$ is approximated by

[TABLE]

with $C_{A}^{t+1}$ the adjacency matrix containing only edges between $\partial\mathcal{G}_{A}^{t+1}$ and $\mathcal{G}_{A}^{t}$. This procedure corresponds to a sum over all paths that are made by moving within $\partial\mathcal{G}_{A}^{t+1}$, shifting back in time to $\mathcal{G}_{A}^{t}$, and then taking any precalculated path in $\mathcal{G}_{A}^{t}$. $C_{A}^{t+1}$ and ${(\mathbf{1}-\partial A_{A}^{t+1})}^{-1}$ can be precalculated, and so the runtime computation requirement is reduced to the product in Eq. (23). This in turn can be sparsified, as $C_{A}^{t+1}$ only contains connections to vertices in $\mathcal{G}_{A}^{t}$ close to the time boundary, and we can delete all terms in $A_{S}^{t}$ that do not connect from these vertices to errors.

We have used the second method for our MWPM decoder, as we expect the error from neglecting higher-order combinations of errors to be small. In order to check this assumption, in Fig. 12 we repeat our simulation protocol with a modified physical error model that excludes all $Y$ and measurement errors. We see that in the absence of these errors, the MWPM decoder performs within the error margin of the decoder upper bound. Note that a small deviation is expected from the discrepancy between a MWPM decoder and a maximum-likelihood decoder Heim16 . With the parameters used in this work, we do not observe any loss of fidelity when we stop accounting for the difference in error rates between ancilla states. We account this to the large error contribution from photon noise and gate infidelity on the ancilla qubits, which do not have this asymmetry. We further note that we operate in a regime of large ancilla error; as described in the text this makes the system counter-intuitively less sensitive to ancilla noise. In systems where this is not the case, it could be that accounting for ancilla asymmetry provides a useful computational method to improve $\epsilon_{\mathrm{L}}$.

Appendix G Implementation of a look-up table decoder

In Tomita14 , the authors describe a decoding scheme specific to Surface-$17$, which is optimized to be implementable with limited computational resources in a short cycle time. This decoding scheme works by using a short decision tree to connect errors to each other in a style similar to blossom. Indeed, this scheme is equivalent to a blossom decoder with all horizontal, vertical and diagonal weights equal Tomita14 . As such, we have implemented the new weights in the blossom decoder rather than utilizing the exact method given.

Appendix H Details of lowest-order approximation

We detail the approximation made to study Surface-$49$ in Sec. II C. Note that this calculation is only for $X$ errors, which are measured by the $Z$ ancillas. This implies that our approximation should attempt to realize the result of blossom, rather than the decoder upper bound.

We begin with the $\mathcal{G}_{A}$ graph defined in App. F. In the absence of correlated errors that cause more than two error signals, any experiment can be approximately described by choosing a set $S\subset E_{A}$ of edges on the graph and assuming the errors that correspond to these edges have occurred. Each ancilla measurement corresponds to a vertex in $\mathcal{G}_{A}$, which records an error if an odd number of edges in $S$ point to the vertex. Each combination $M_{a}$ of ancilla measurements can be generated by multiple error sets $S$.

Formally, let us write $\mathcal{M}$ for the set of all combinations of ancilla measurements and $\mathcal{S}$ for the set of all combinations of errors (so $\mathcal{S}=2^{E_{A}}$). We then define a function $\phi:\mathcal{S}\rightarrow\mathcal{M}$ that takes a combination of errors to the resultant measurement outcomes. Let us fix a logical $Z$ operator $Z_{L}$ on the surface-code fabric. Then to each $S\in\mathcal{S}$ we can assign a parity $p(S)=\pm 1$ depending on whether the product of all errors in $S$ commute with $Z$ or not. A decoding then consists of a choice of parity $p_{d}(M)$ for each $M\in\mathcal{M}$. Such a decoding correctly decodes $S\in\mathcal{S}$ if $p_{d}(\phi(S))=p(S)$, and creates a logical error otherwise. The source of logical errors in a perfect decoder is then precisely the fact that we can have two error combinations $S_{1},S_{2}\in\mathcal{S}$ such that $\phi(S_{1})=\phi(S_{2})$ but $p(S_{1})\neq p(S_{2})$.

The above suggests a method by which a perfect decoder can be constructed. As defined, $\phi^{-1}(M)\subset\mathcal{S}$ is the set of error combinations $S$ that return a measurement $M_{\in}\mathcal{M}$. For each error combination $S$, we can calculate the probability of this occurring:

[TABLE]

The optimal choice of $p_{d}(M_{a})$ is the one maximizing

[TABLE]

and the fidelity of such a decoder (over the entire experiment) can be calculated as

[TABLE]

At this point the only approximation that has been made is to neglect the $T_{1}$ asymmetry in the system, which we have shown previously in this work to be negligible. Unfortunately, the above function cannot be evaluated exactly; the number of error combinations $S$ is approximately $2^{200}$ for $4$ cycles of Surface-$49$. Our goal instead is to approximate this to the lowest order in the physical qubit error rate.

Let us make the approximation that our error combinations $S$ can be split into small, well-separated pieces of errors containing separate correctable and non-correctable parts, $S=\cup_{i}S^{i}$. To each $S^{i}$ we can assign a time step $t(S^{i})$, being the earliest time of the first error measurement observed (in $\phi(S^{i})$). The error rate per round, $\epsilon_{\mathrm{L}}$, can be determined by summing Eq. 26 over all pieces $S^{i}$ of all combinations $S$ such that $t(S^{i})=T$ (with arbitrary $T$), as the effect of repeated errors from $S^{i},S^{j}\subset S$ is taken into account during the derivation of the logical fidelity equation (Eq. 2 in the main text).

[TABLE]

Let us also extend the above division of $S$ to a division of $M$ into separate pieces $M^{a}$, and rewrite our sum slightly,

[TABLE]

Where here we have brought the sum over the global combinations of syndromes and measurements inside a new function $\bar{r}$

[TABLE]

If we took this approximation literally and considered the sum over every possible combination $S$ containing $S^{i}$, the final sum in Eq. 29 would reduce to

[TABLE]

However, this includes error combinations $S$ that cannot be easily separated into $S^{i}$ and ‘something else’, i.e. they contain other errors $e$ that cannot be separated from $S^{i}$. Eq. 30 is then equivalent to assuming that if $S^{i}$ is an uncorrectable logical error, no nearby combination of physical errors $S^{\prime}$ can be combined such that $S^{i}\cup S^{\prime}$ is correctable unless $S^{\prime}$ itself is an uncorrectable logical error. Such combinations would serve to reduce the calculated $\epsilon_{L}$, and so $\bar{r}^{(u)}$ gives an upper bound for $\epsilon_{L}$ in Eq. 28. For a lower bound, we approximate that for any uncorrectable error combination $S^{i}$, approximately one rounds-worth of single errors would undo the logical error, leading to the approximation

[TABLE]

We now make one further approximation, and sum Eq. 28 only over the shortest $S^{i}$ that can be expected to contribute to the final error rate. That is, we sum over those $S^{i}$ with $|S^{i}|\leq(d+1)/2$, and that spread directly across the chain. The error incurred from this approximation is roughly proportional to the largest single error, which is no more than $5\%$ throughout our study. We use $\bar{r}^{(u)}$ and $\bar{r}^{(l)}$ to give the error bars shown in Fig. 5. Points in the plot are taken as a log average of the upper and lower bounds, and thus have no particular relevance themselves. We see that the numerical calculation falls within the corresponding error bars for almost the entire dataset, giving verification for our method, save a slight deviation at one point where it falls below. Moreover, as the simulated Surface-$17$ error rate lies above the upper bound found for the Surface-$49$ error rate (with the standard set of parameters from the main text), our claim that Surface-$17$ will operate below the fault-tolerant threshold is quite strong.