Efficient simulation of quantum error correction under coherent error based on non-unitary free-fermionic formalism

TL;DR

This paper introduces an efficient method using non-unitary free-fermionic formalism to accurately simulate quantum error correction under coherent noise, revealing the impact of noise coherence on error thresholds.

Contribution

It develops a novel, efficient simulation scheme for quantum error correction under non-Clifford, coherent noise, enabling precise threshold evaluation without approximation.

Findings

Error threshold drops to one third with fully coherent noise.

Method applicable to 1D repetition and surface codes.

Threshold dependence on coherence explained by leading-order analysis.

Abstract

In order to realize fault-tolerant quantum computation, tight evaluation of error threshold under practical noise models is essential. While non-Clifford noise is ubiquitous in experiments, the error threshold under non-Clifford noise cannot be efficiently treated with known approaches. We construct an efficient scheme for estimating the error threshold of one-dimensional quantum repetition code under non-Clifford noise. To this end, we employ non-unitary free-fermionic formalism for efficient simulation of the one-dimensional repetition code under coherent noise. This allows us to evaluate the effect of coherence in noise on the error threshold without any approximation. The result shows that the error threshold becomes one third when noise is fully coherent. Our scheme is also applicable to the surface code undergoing a specific coherent noise model. The dependence of the error threshold on noise coherence can be explained with a leading-order analysis with respect to coherence terms in the noise map. We expect that this analysis is also valid for the surface code since it is a two-dimensional extension of the one-dimensional repetition code. Moreover, since the obtained threshold is accurate, our results can be used as a benchmark for approximation or heuristic schemes for non-Clifford noise.

Figures

Equations

$$\mathcal{E}(\rho ,X)$$

$$\mathrm{Pr}({\mathbf{t}}_k)=\Gamma ({\mathbf{t}}_k):=⟨0^{\otimes n}⟩({\mathbf{K}}^{({\mathbf{t}}_k)}{)}^{†}{\mathbf{K}}^{({\mathbf{t}}_k)}{0}^{\otimes n},$$

$$\mathrm{Pr}(a_{\mathrm{out}}$$

$$p_{\mathrm{L}}=\mathrm{Pr}(a_{\mathrm{out}}=1)={⟨{\Gamma }_{\mathrm{L}}(\mathbf{t})/\Gamma (\mathbf{t})⟩}_{\mathbf{t}},$$

$$\Gamma ({\mathbf{t}}_k)=⟨\tilde{\psi }⟩({\mathbf{K}}^{({\mathbf{t}}_k)}{)}^{†}{\mathbf{K}}^{({\mathbf{t}}_k)}\tilde{\psi },$$

$${\Gamma }_{\mathrm{L}}(\mathbf{t})=⟨\tilde{\psi }⟩(R(\mathbf{t}){\mathbf{K}}^{(\mathbf{t})}{)}^{†}\frac{I-Z_n{Z}_{n+1}}{2}R(\mathbf{t}){\mathbf{K}}^{(\mathbf{t})}\tilde{\psi },$$

$$\tilde{M}$$

$$M({\mathbf{t}}_k)$$

$$\Gamma ({\mathbf{t}}_k)$$

$$\frac{{\Gamma }_{\mathrm{L}}(\mathbf{t})}{\Gamma (\mathbf{t})}$$

$$p_{\mathrm{L}}(d,p)\sim \alpha {\left(\frac{p}{p_{\mathrm{th}}}\right)}^{d/2},$$

$$\lambda (d,p):=\frac{p_{\mathrm{L}}(d+2,p)}{p_{\mathrm{L}}(d,p)}\sim \frac{p}{p_{\mathrm{th}}},$$

$$\rho (\mathcal{C}):=(\mathcal{C}\otimes \mathrm{Id})(|{\varphi }_{\mathrm{init}}⟩⟨{\varphi }_{\mathrm{init}}|)$$

$$|{\varphi }_{\mathrm{init}}⟩:=\frac{|0_{\mathrm{L}}⟩|0⟩+|1_{\mathrm{L}}⟩|1⟩}{\sqrt{2}},$$

$$c_{2i-1}$$

$$c_{2i}$$

$${\tilde{S}}_i^{\mathrm{r},c}=(\prod _{j=i}^{d-1}{S}_j^{\mathrm{r},c})S^{\mathrm{g},c}$$

$$W_{d(c-1)+j}{X}_{d(c-1)+j+1}\cdots X_{d(c+1)-j}{W}_{d(c+1)+1-j},$$

$$S_i^{\mathrm{r},c}=\{\begin{array}{ll}{\tilde{S}}_i^{\mathrm{r},c}{\tilde{S}}_{i+1}^{\mathrm{r},c}& \text{if }i<d-1\\ {\tilde{S}}_{d-1}^{\mathrm{r},c}{S}^{\mathrm{g},c}& \text{if }i=d-1\end{array}.$$

$$P(S,s):=\frac{1}{2}(I+(-1{)}^{s}S).$$

$$P(S^{\mathrm{g},c},s^{\mathrm{g}})\prod _{i=1}^{d-1}P(S_i^{\mathrm{g},c},s_i^{\mathrm{r}})=P(S^{\mathrm{g},c},s^{\mathrm{g}})\prod _{i=1}^{d-1}P({\tilde{S}}_i^{\mathrm{g},c},{\tilde{s}}_i^{\mathrm{r}})$$

$$s_i=\{\begin{array}{ll}({\tilde{s}}_i^{\mathrm{r}}+{\tilde{s}}_{i+1}^{\mathrm{r}})\mathrm{mod}2& \text{if }j<d-1\\ ({\tilde{s}}_{d-1}^{\mathrm{r}}+s^{\mathrm{g}})\mathrm{mod}2& \text{if }j=d-1\end{array}.$$

$${\mathcal{P}}_{\mathrm{C}}$$

$$L_Z$$

$$L_Y$$

$$|{\varphi }_{\mathrm{init}}⟩⟨{\varphi }_{\mathrm{init}}|={\mathcal{P}}_{\mathrm{C}}P(L_Z{Z}_{n+1},0)P(-L_Y{Y}_{n+1},0).$$

$$p_{s^{\prime }}:=⟨\varphi ⟩P(S,s^{\prime })\varphi (s^{\prime }=0,1),$$

$$\frac{P(S,s^{\prime })|\varphi ⟩}{\sqrt{⟨\varphi ⟩P(S,s^{\prime })\varphi }}.$$

$$L_X=-i{L}_Y{L}_Z=-Y_d\left(\prod _{i=d+1}^{n-d}{X}_i\right)Y_{n-d+1}.$$

$$|{\varphi }_{\mathrm{final}}⟩⟨{\varphi }_{\mathrm{final}}|$$

$$A_{W,W^{\prime }}$$

$$\rho (\mathcal{C})=\mathrm{ave}(Z_{n+1}^w|{\varphi }_{\mathrm{final}}⟩⟨{\varphi }_{\mathrm{final}}|Z_{n+1}^w),$$

$$F(\mathcal{C})$$

Full text

Introduction.— Quantum error correction (QEC) is a key technology for building a scalable fault-tolerant quantum computer. According to the theory of fault-tolerant quantum computation, one can perform quantum computation with arbitrary accuracy if the error probability is below a certain threshold value Kitaev (1997); Knill et al. (1998); Aharonov and Ben-Or (1997). The threshold values of various QEC schemes have been calculated under various assumptions of the noise models and degrees of rigor Fern et al. (2006); Greenbaum and Dutton (2016); Wang et al. (2003, 2011); Fowler et al. (2012a); Stephens (2014); Ghosh et al. (2012); Geller and Zhou (2013); Gutiérrez et al. (2013); Gutiérrez and Brown (2015); Magesan et al. (2013); Puzzuoli et al. (2014); Gutiérrez et al. (2016); Rahn et al. (2002); Chamberland et al. (2016); Tomita and Svore (2014); Ferris and Poulin (2014); Darmawan and Poulin (2016). In the case of the noise model which only consists of probabilistic Clifford gates and Pauli measurement channels, such as depolarizing noise, the threshold value can be efficiently and accurately estimated numerically Wang et al. (2003, 2011); Fowler et al. (2012a); Stephens (2014) by virtue of the Gottesman-Knill theorem Gottesman (1998); Aaronson and Gottesman (2004). On the other hand, non-Clifford noise is unavoidable in practical experiments Kelly et al. (2015); Córcoles et al. (2015); Ristè et al. (2015), but QEC circuits under non-Clifford noise cannot be treated with this approach. Specifically, it is theoretically predicted that coherent noise, which is non-Clifford and is caused, for example, by over rotation, can have negative effects on quantum error correction Kueng et al. (2016). Therefore, massive effort has been made for evaluating the effect of noise coherence on the error threshold. Since the simulation of quantum circuits under arbitrary local noise sometimes becomes as hard as that of universal quantum computation, we cannot efficiently simulate QEC circuits under coherent noise with straightforward methods. While its computational cost can be relaxed in some extent Tomita and Svore (2014); Ferris and Poulin (2014); Darmawan and Poulin (2016), the tractable number of qubits with straightforward methods is limited. In the case of concatenated codes, there is an efficient method to analytically estimate the error threshold under non-Clifford noise Rahn et al. (2002). However, this technique is not applicable to topological codes, which are more feasible in practical experiments Kelly et al. (2015); Córcoles et al. (2015); Ristè et al. (2015); Lidar and Brun (2013). In general, we may approximate non-Clifford noise by a Clifford channel for an efficient simulation Gutiérrez et al. (2016); Ghosh et al. (2012); Geller and Zhou (2013); Gutiérrez et al. (2013); Gutiérrez and Brown (2015); Magesan et al. (2013); Puzzuoli et al. (2014), but the accuracy of the estimated threshold is sacrificed. An efficient and accurate scheme, which works for topological codes under non-Clifford noise, is still lacking.

Here, we construct an efficient and accurate scheme to simulate one-dimensional (1D) repetition code with repetitive parity measurements under coherent noise. While the 1D repetition code cannot protect a logical qubit from arbitrary single-qubit error, it is still able to capture a necessary ingredient for fault-tolerant QEC, and hence was experimentally demonstrated as a building block for scalable fault-tolerant quantum computation Kelly et al. (2015); Ristè et al. (2015). The key idea in our scheme is reducing the QEC circuit of the 1D repetition code to a classically simulatable class of non-unitary free-fermionic dynamics, which is known as a variation of matchgate quantum computing Valiant (2002); Terhal and DiVincenzo (2002); Knill (2001); Bravyi (2005a, b); Jozsa and Miyake (2008); Jozsa et al. (2015); Brod (2016); Bravyi et al. (2014). As compared to the stochastic noise model, we find that the error threshold of the 1D repetition code becomes about one third when noise is fully coherent. The dependence of the error threshold on noise coherence is explained by using a leading-order analysis with respect to coherence terms of noise map. We expect that a similar analysis holds in the surface code Dennis et al. (2002); Bravyi and Kitaev (1998); Fowler et al. (2012b), which is the most experimentally feasible QEC scheme, since the surface code is a two-dimensional extension of the 1D repetition code.

Furthermore, our accurate results can be used as a benchmark for approximation or heuristic schemes for estimating the error threshold under non-Clifford noise.

Simulation of the 1D repetition code.— The quantum circuit of the 1D repetition code with repetitive parity measurements is shown in Fig. 1(a). In the 1D repetition code, one logical bit $a_{\rm in}$ is encoded into physical $n$ data qubits $\{\ket{0}^{\otimes n},\ket{1}^{\otimes n}\}$, which are stabilized by operators $\{Z_{i}Z_{i+1}\}_{i=1}^{n-1}$, where $A_{i}(A\in\{X,Y,Z\})$ is the Pauli operator on the $i$-th data qubit. Error syndrome is measured through $(n-1)$ measurement qubits, each of which monitors the parities of the neighboring data qubits, i.e., $Z_{i}Z_{i+1}$. The measurements are repetitively performed for $T$ cycles. The encoded bit is finally decoded from the data qubits and $T(n-1)$ syndromes, which can be efficiently done using minimum-weight perfect matching Kelly et al. (2015). The probability with which the decoded bit $a_{\rm out}$ is flipped is defined as logical error probability $p_{\rm L}$, which is the failure probability of the decoding.

Since the 1D repetition code is capable of correcting only $X$-type error, we consider a CPTP (completely positive trace-preserving) map of a general single-qubit $X$-type noise, which is regarded as a mixture of the $X$-type unitary (fully-coherent) and stochastic (incoherent) noise:

[TABLE]

where $\theta$ is defined by $\cos\theta=\sqrt{1-p}$ and $\sin\theta=\sqrt{p}$. The parameter $c$ ($0\leq c\leq 1$), which we call noise coherence, is a measure of coherence in the noise. We call the parameter $p$ ($0\leq p\leq 1$) as the physical error probability since it can be understood as the probability with which the input state $\ket{0}$ is measured as the output state $\ket{1}$. We consider two types of noise allocation models Landahl et al. (2011) as shown in Fig. 1(b). In the case of the phenomenological model, a noise map ${\mathcal{E}(\rho,X)}$ is located on each of the data and measurement qubits at the beginning of each cycle. In the case of the circuit-based model, the noise map is located at each time step of preparation, gate operation, and measurement, on every qubit including the one that is idle at the time step. There, we assume that two-qubit noise map $p_{XI}{\mathcal{E}}(\rho,X\otimes I)+p_{IX}{\mathcal{E}}(\rho,I\otimes X)+p_{XX}{\mathcal{E}}(\rho,X\otimes X)$ acts on the output qubits after each controlled-Not (CNOT) operation. Since our noise models are symmetric over the bit values, we may choose $a_{\rm in}=0$ to evaluate the logical error probability as $p_{\rm L}={\rm Pr}(a_{\rm out}=1)$.

Before reducing the noisy circuit to free-fermionic dynamics, we reformulate it as a sequence of generalized measurements on the data qubits, such that the state after each measurement is pure. We denote the outcome of the $k$-th measurement by $t_{k}$ and the corresponding Kraus operator by $K_{k}^{(t_{k})}$. The probability of a sequence of outcomes ${\bm{t}}_{k}:=t_{k}...t_{1}$ is given by

[TABLE]

where ${\bm{K}}^{({\bm{t}}_{k})}:=K_{k}^{(t_{k})}K_{k-1}^{(t_{k-1})}...K_{1}^{(t_{1})}$. We may identify three types of operations on data qubits to assign Kraus operators. For clarity, we describe the case of the phenomenological model. The first type is the single-qubit noise ${\mathcal{E}}$ given in Eq.(1). Its operation on the $i$-th qubits is equivalently described by Kraus operators $K_{{\rm noise},i}^{(\phi)}=\sqrt{p(\phi)}e^{i\phi X_{i}}$, where $\phi\in\{\pm\theta\}$ and $p(\pm\theta):=\frac{1\pm c}{2}$. The second type is the parity measurement on the $i$-th and $(i+1)$-th data qubits, which composed of a measurement qubit and two CNOT gates. Treating the noise map $\mathcal{E}$ on the measurement qubit as above, it can be represented by $K^{(s,\phi)}_{{\rm parity},i}=\sqrt{p(\phi)}\frac{1}{2}(I+(-1)^{s}e^{-2i\phi}Z_{i}Z_{i+1})$, where $s\in\{0,1\}$ is the output of the parity measurement. In the case of the circuit-based model, we may still use the same form of $K^{(s,\phi)}_{{\rm parity},i}$ except for varying the probability mass function $p(\phi)$ (see Appendix A SM ). The third type appears in an alternate description of the decoding process. Though the input bit is usually decoded through noisy direct measurements of the data qubits and classical computation, we use the following equivalent process instead. We apply map $\mathcal{E}$ on each data qubit. We perform ideal parity measurements on neighboring data qubits, whose Kraus operator is given by $K^{(s)}_{{\rm parity},i}=\frac{1}{2}(I+(-1)^{s}Z_{i}Z_{i+1})$. Let $k_{f}$ be the index of the last ideal parity measurement, and define ${\bm{t}}:={\bm{t}}_{k_{f}}$. Based on all the measured parities, which is included in $\bm{t}$, we choose a recovery operation $R({\bm{t}})=\prod_{i=1}^{n}X_{i}^{r_{i}({\bm{t}})}$, where $r_{i}({\bm{t}})\in\{0,1\}$ is determined using minimum-weight perfect matching Kelly et al. (2015); Edmonds (1965); Kolmogorov (2009). The recovered state $R({\bm{t}}){\bm{K}}_{k_{f}}^{({\bm{t}})}\ket{0}^{\otimes n}$ is in the code space of the 1D repetition code, which should be written in the form $\alpha\ket{0}^{\otimes n}+\beta\ket{1}^{\otimes n}$. The decoded bit $a_{\rm out}$ is thus obtained by measuring the $n$-th qubit. The joint probability of obtaining $\bm{t}$ and a decoding failure is then given by

[TABLE]

From Eqs. (2) and (3), we have,

[TABLE]

which means that we can accurately calculate $p_{\rm L}$ by sampling $\bm{t}$ with probability $\Gamma(\bm{t})$ repeatedly and by taking the average of $\Gamma_{\rm L}({\bm{t}})/\Gamma({\bm{t}})$. Since the sampling of $\bm{t}$ can be done by sequentially generating $t_{k}$ according to ${\rm Pr}(t_{k}|{\bm{t}}_{k-1})=\Gamma({\bm{t}}_{k})/\Gamma({\bm{t}}_{k-1})$, the efficiency of this scheme follows that of computing $\Gamma({\bm{t}}_{k})$ and $\Gamma_{\rm L}(\bm{t})$.

Reduction to non-unitary free-fermionic dynamics.— We use non-unitary free-fermionic dynamics to calculate $\Gamma({\bm{t}}_{k})$ and $\Gamma_{\rm L}({\bm{t}})$ efficiently. Let us briefly summarize the known facts about non-unitary free-fermionic dynamics Bravyi (2005a); Bravyi et al. (2014). We define $\{c_{i}\}$ ($1\leq i\leq 2n$) as the Majorana fermionic operators for $n$ fermionic modes, which satisfy $\{c_{i},c_{j}\}=2\delta_{i,j}$, $c_{i}^{\dagger}=c_{i}$ and $c_{i}^{2}=I$. The covariance matrix for a pure state $\ket{\psi}$ is defined as $M_{ij}=\frac{-i}{\braket{\psi}{\psi}}\frac{\braket{\psi}{[c_{i},c_{j}]}{\psi}}{2}$. We call the state $\ket{\psi}$ is a fermionic Gaussian state (FGS) iff the covariance matrix satisfies $MM^{T}=I$. An FGS can be fully specified by a pair $(M,\Gamma)$, where $\Gamma$ is the norm $\braket{\psi}{\psi}$. The absolute value of the inner product of two FGSs, $\ket{\psi}\mapsto(M_{\psi},\Gamma_{\psi})$ and $\ket{\phi}\mapsto(M_{\phi},\Gamma_{\phi})$, is given by $|\braket{\psi}{\phi}|^{2}=2^{-n}\Gamma_{\psi}\Gamma_{\phi}\det(M_{\psi}+M_{\phi})$. An operator of form $e^{\sum_{i<j}\alpha_{ij}c_{i}c_{j}}$ with $\alpha_{ij}$ being a complex value is called a fermionic Gaussian operator (FGO). Note that FGOs are not necessarily unitary. An FGO maps any FGS to another FGS. Given an FGO $G$ and an input FGS $\ket{\psi}\mapsto(M,\Gamma)$, the description $(M^{\prime},\Gamma^{\prime})$ for the output state $G\ket{\psi}$ is calculated as follows. Consider a fermionic maximally entangled state $\ket{\psi_{\rm M}}\mapsto(M_{\rm M},\Gamma_{\rm M})$ of $2n$ fermionic modes, which is defined by M_{\rm M}=\left(\begin{array}[]{cc}0&I\\ -I&0\end{array}\right) and $\Gamma_{\rm M}=1$. We calculate $(M_{\rm G},\Gamma_{\rm G})$ corresponding to the state $\ket{\psi_{\rm G}}=(G\otimes I)\ket{\psi_{\rm M}}$. In terms of matrices $A,B,D$, which are defined by \left(\begin{array}[]{cc}A&B\\ -B^{T}&D\end{array}\right)=M_{\rm G}, the output is calculated as $(M^{\prime},\Gamma^{\prime})=(A-B(M-D)^{-1}B^{T},\Gamma_{\rm G}\Gamma\sqrt{\det(M-D)})$. In this way, free-fermionic dynamics consisting of FGOs on FGSs can be simulated efficiently.

Now we are ready to reformulate the QEC process with non-unitary free-fermionic dynamics. Using the Jordan-Wigner transformation, we may choose $c_{2i-1}=\left(\prod_{j=1}^{i-1}X_{j}\right)Z_{i}$ and $c_{2i}=\left(\prod_{j=1}^{i-1}X_{j}\right)Y_{i}$. We see that $X_{i}=-ic_{2i}c_{2i-1}$ and $Z_{i}Z_{i+1}=-ic_{2i+1}c_{2i}$, which are both quadratic terms of the Majorana fermionic operators. Therefore, $K_{{\rm noise},i}^{(\phi)},K_{{\rm parity},i}^{(s,\phi)},K_{{\rm parity},i}^{(s)}$ and $R$ are FGOs. Unfortunately, the initial state $\ket{0}^{\otimes n}$ and the operator $\frac{I-Z_{n}}{2}$ for calculating $\Gamma_{\rm L}({\bm{t}})$ are not an FGS and an FGO, respectively. For an efficient simulation, we need a further trick as follows. (While similar tricks in Refs. Bravyi (2005b); Jozsa et al. (2015); Brod (2016) might be employed, the following construction is much simpler and more efficient for our purpose.) We add the $(n+1)$-th ancillary qubit and corresponding Majorana fermionic operators $c_{2n+1},c_{2n+2}$. Using an FGS $\ket{\tilde{\psi}}:=(\ket{0}^{\otimes(n+1)}+\ket{1}^{\otimes(n+1)})/\sqrt{2}$, it is not difficult to show that

[TABLE]

since $K^{({\bm{t}}_{k})}$ and $R({\bm{t}})$ commute with $\prod_{j=1}^{n+1}X_{j}$. Hence, they can be efficiently calculated (see Appendix A SM for detail).

Result.— We show the logical error probability $p_{\rm L}$ as a function of the physical error probability $p$ under incoherent noise ($c=0$) and fully coherent noise ($c=1$) in Fig. 2. To observe clear behavior of the error threshold, we have varied the number $T$ of cycles according to $n$ as $T=n-1$. We also assumed uniform error probability for two-qubit noise, i.e. $p_{XI}=p_{IX}=p_{XX}=\frac{1}{3}$. We employed uniform weighting for performing minimum-weight perfect matching. The logical error probability $p_{\rm L}$ is expected to be exponentially small in the number of the data qubits $n$ as far as the physical error probability $p$ is below a certain value, which we call the error threshold $p_{\rm th}$. By using the scaling ansatz Wang et al. (2003); Stephens (2014), we obtained the threshold values $p_{\rm th}=10.34(1)\%$ for $c=0$ and $7.87(2)\%$ for $c=1$ in the phenomenological model, and $3.243(6)\%$ for $c=0$ and $1.040(5)\%$ for $c=1$ in the circuit-based model. Our result for $c=0$ in the case of the phenomenological model is consistent with the known results Wang et al. (2003). For more detailed procedures, see Appendix B SM . We also confirmed exponential decay of logical error probability $p_{\rm L}$ with code distance $d$ below the threshold value, which is approximated by $p_{\rm L}\propto\left(\frac{p}{p_{\rm th}}\right)^{d/2}$ regardless of the coherence of the noises (see Appendix E SM ).

Dependence of the error threshold $p_{\rm th}$ on the noise coherence $c$ is shown in the insets of Fig. 2. We see that the error threshold $p_{\rm th}$ decreases as the noise coherence $c$ increases. Note that non-uniform weighting improves the error threshold, but only slightly (see Appendix C SM ). The dependence on $c$ can be explained with a leading-order analysis as follows. The noise map of Eq.(1) can be rewritten as ${\mathcal{E}}(\rho,X)=(1-p)\rho+ic\sqrt{(1-p)p}(X\rho-\rho X)+pX\rho X$. We call the second term $ic\sqrt{(1-p)p}(X\rho-\rho X)$ as the coherence term. This term contributes to diagonal terms of the density matrix only through a concatenation of multiple noise maps. The correction to the diagonal terms after several cycles is written as even-order terms in $c\sqrt{(1-p)p}$. For $p\ll 1$, the leading order of the correction is $O(p)$ in the circuit-based model, while it is $O(p^{2})$ in the phenomenological model since an error on a data qubit spreads to two measurement qubits before the next noise map is applied on the data qubit. For example, the product of the coherence terms of noise maps located in the positions shown in Fig. 3 contributes the correction. In the case of the phenomenological model, the leading term is proportional to $c^{4}p^{2}$, and its sign depends on the results of previous syndrome measurements. Such a noise leads to space-time correlations in the syndrome measurements. Since the decoder is not adapted to such correlations, the existence of coherence in noise is expected to result in a worse logical error probability. On the other hand, in the case of the circuit-based model, the leading term is proportional to $c^{2}p$, and it always increases the error probability. This directly worsens the logical error probability and the error threshold.

In the case of the circuit-based model, we seek a more quantitative explanation of the behavior by proposing a heuristic ansatz as follows. We define an effective physical error probability of a data qubit per cycle $p_{\rm eff}(p,c)$ (the precise definition is given in Appendix D SM ). The probability $p_{\rm eff}(p,c)$ should be expanded for small $p$ as $p_{\rm eff}(p,c)=\beta(1+\alpha c^{2})p+O(p^{2})$, where $\alpha$ is constant and is independent of the system size $n$. We assume that the logical error probability $p_{\rm L}(p,c)$ can be well explained by the local increase of noise, i.e., $p_{\rm L}(p,c)=p_{\rm L}((1+\alpha c^{2})p,0)$. Based on this ansatz, the error threshold under coherent noise $p_{\rm th}(c)$ can be written as $p_{\rm th}(c)\sim\frac{p_{\rm th}(0)}{1+\alpha c^{2}}$ if $p_{\rm th}\ll 1$. By using the analytically obtained value $\alpha=11/6$ (see Appendix D SM ), this ansatz gives the solid curve in the inset of Fig. 2, which is in good agreement with the accurate numerical results. We may expect that a similar leading-order ansatz also holds for the surface code, since it is a two-dimensional extension of the 1D repetition code. The factor $\alpha$ is also easily obtained by analytically calculating the effective bit-flip probability, and $p_{\rm th}(0)$ for incoherent noises can be efficiently computed. Therefore, the error threshold of the surface code under coherent noise will be estimated by the same approach.

Conclusion and discussion.— We constructed an efficient and accurate scheme for estimating the error threshold of the 1D repetition code under coherent noise. We have calculated the error threshold under coherent noise in terms of the physical error probability $p$ and the noise coherence $c$. The parameters $p$ and $c$ can be experimentally accessible by randomized and purity benchmarkings, respectively Knill et al. (2008); Wallman et al. (2015). We emphasize here that the proposed accurate and efficient scheme is not limited to the 1D repetition code. In fact, in Appendix F SM we provide an example with a fully quantum code, which simulates the surface code under a phenomenological coherent noise model. We have also proposed a leading-order ansatz for the estimation of the error threshold under coherent noise, and found that it reproduces the accurate numerical results well. This suggests that the effect of the coherent noise on the surface code will be assessed by an analogous ansatz, which can be calculated easily. In more general terms, the obtained accurate error thresholds of the 1D repetition code will serve as a reference to test the accuracy of approximation or heuristic schemes for simulating non-Clifford noise, as was done for the leading-order ansatz.

Efficient simulation of quantum error correction under coherent error based on non-unitary free-fermionic formalism - Supplemental material

Appendix A Appendix A: Detail of sampling scheme

In this appendix, we describe the scheme of sampling $\bm{t}$ and computing $\Gamma_{\rm L}(t)/\Gamma(t)$ in detail. The simulation can be divided into three processes.

*Process 1 — Allocation of single-qubit noise maps

In the case of the phenomenological model, the allocation of the noise maps is fixed and follows Fig. 1(b). In the case of the circuit-based model, it is probabilistically chosen for each sampling as follows. The two-qubit noise map assumed after each CNOT gate contains a non-local $X$-type noise term $e^{iX_{\rm control}X_{\rm target}}$, which is not an FGO. We can convert this non-local noise to local noise by replacing it with a local noise preceding the CNOT gate $U_{\rm CNOT}$ since ${\mathcal{E}}(U_{\rm CNOT}\rho U_{\rm CNOT}^{\dagger},X_{\rm control}X_{\rm target})=U_{\rm CNOT}{\mathcal{E}}(\rho,X_{\rm control}I_{\rm target})U_{\rm CNOT}^{\dagger}$. Thus, to simulate the two-qubit noise map faithfully, we probabilistically place a single-qubit noise map ${\mathcal{E}}$ at 1) the target qubit after the CNOT gate, 2) the control qubit after the CNOT gate, and 3) the control qubit before the CNOT gate with the probabilities $p_{XI},p_{IX}$ and $p_{XX}$, respectively. The single-qubit noises associated with state preparations and measurements for the measurement qubits are placed deterministically. The ones for the data qubits at the beginning of the decoding process are also placed deterministically.

*Process 2 — Simulation of the circuit

The covariance matrix of the state $\ket{\tilde{\psi}}=(\ket{0}^{\otimes(n+1)}+\ket{1}^{\otimes(n+1)})/\sqrt{2}$ is given as follows.

[TABLE]

We start the simulation from $(\tilde{M},1)$, which is formally denoted by $(M({\bm{t}}_{0}),\Gamma({\bm{t}}_{0}))$. Given $(M({\bm{t}}_{k-1}),\Gamma({\bm{t}}_{k-1}))$, the value of $t_{k}$ is sampled from the probability $\Gamma({\bm{t}}_{k})/\Gamma({\bm{t}}_{k-1})$. Then the updated pair $(M({\bm{t}}_{k}),\Gamma({\bm{t}}_{k}))$ is calculated by

[TABLE]

where $A^{(t_{k})},B^{(t_{k})},D^{(t_{k})}$ and $\Gamma_{\rm G}^{(t_{k})}$ are associated with the FGO $K_{k}^{(t_{k})}$.

There are three types of operators $K_{{\rm noise},i}^{(\phi)},K_{{\rm parity},i}^{(s,\phi)}$ and $K_{{\rm parity},i}^{(s)}$ for $K_{k}^{(t_{k})}$. For $K_{{\rm noise},i}^{(\phi)}=\sqrt{p(\phi)}e^{i\phi X_{i}}$ with $p(\pm\theta)=\frac{1\pm c}{2}$, which is a noise operation on a data qubit, we have $\Gamma_{\rm G}^{(t_{k})}=p(\phi)$ and $A^{(t_{k})}=D^{(t_{k})}=0$. The submatrix $B^{(t_{k})}$ is calculated as $B^{(t_{k})}=I+B^{\prime}$, where $B^{\prime}$ has nonzero elements only for \left(\begin{array}[]{cc}B^{\prime}_{2i-1,2i-1}&B^{\prime}_{2i-1,2i}\\ B^{\prime}_{2i,2i-1}&B^{\prime}_{2i,2i}\end{array}\right)=\left(\begin{array}[]{cc}-1+\cos 2\phi&-\sin 2\phi\\ \sin 2\phi&-1+\cos 2\phi\end{array}\right).

For $K_{{\rm parity},i}^{(s,\phi)}=\sqrt{p(\phi)}\frac{1}{2}(I+(-1)^{s}e^{-2i\phi}Z_{i}Z_{i+1})$, which represents a parity measurement on two qubits, we have $\Gamma_{\rm G}^{(t_{k})}=p(\phi)/2$. The matrix $A^{(t_{k})}$ has nonzero elements only for $A_{2i,2i+1}=-A_{2i+1,2i}=-(-1)^{s}\cos 2\phi$, and $D^{(t_{k})}=-A^{(t_{k})}$. The matrix $B^{(t_{k})}$ is calculated as $B^{(t_{k})}=I+B^{\prime\prime}$, where $B^{\prime\prime}$ has nonzero elements only for \left(\begin{array}[]{cc}B^{\prime\prime}_{2i,2i}&B^{\prime\prime}_{2i,2i+1}\\ B^{\prime\prime}_{2i+1,2i}&B^{\prime\prime}_{2i+1,2i+1}\end{array}\right)=\left(\begin{array}[]{cc}-1&(-1)^{s}\sin 2\phi\\ -(-1)^{s}\sin 2\phi&-1\end{array}\right). The form of $p(\phi)$ depends on the noise allocation model. In the case of the phenomenological model, there is one $X$-type noise map on each measurement qubit, and the probability mass function $p(\phi)$ is defined as $p(\pm\theta)=\frac{1\pm c}{2}$. In the case of the circuit-based model, there are multiple $X$-type noise maps on the measurement qubit. We denote the number of the noise maps as $N$, which depends on the allocation at Process 1. Since a CNOT gate commutes with $X$-type noise map on the target qubit, we are allowed to use the same form of $K_{{\rm parity},i}^{(\phi)}$ except that the probability mass function $p(\phi)$ is replaced by $p(k\theta)=\binom{N}{(N-k)/2}\left(\frac{1+c}{2}\right)^{(N-k)/2}\left(\frac{1-c}{2}\right)^{(N+k)/2}$ for $k=\{-N,-N+2,...,N-2,N\}$.

Finally, for $K_{{\rm parity},i}^{(s)}$, we have $\Gamma_{\rm G}^{(t_{k})}=\frac{1}{2}$, and $A^{(t_{k})},B^{(t_{k})}$ and $D^{(t_{k})}$ are equivalent to those for the $K_{{\rm parity},i}^{(s,\phi)}$ with $\phi=0$.

After repeating the $(n-1)T$ syndrome measurements, we perform noiseless parity measurement for each neighboring data qubits. As a result, we obtain $(n-1)(T+1)$ outputs of the parity measurements $\bm{s}$ included in $\bm{t}$, and the final state of data qubits $(M({\bm{t}}),\Gamma({\bm{t}}))=(M({\bm{t}}_{k_{f}}),\Gamma({\bm{t}}_{k_{f}}))$.

*Process 3 — Decoding

We determine the recovery operation $R({\bm{t}})=\prod_{i=1}^{n}X_{i}^{r_{i}({\bm{t}})}$ by using minimum-weight perfect matching. We then calculate $\Gamma_{\rm L}({\bm{t}})/\Gamma({\bm{t}})$ from Eqs. (5) and (6) as

[TABLE]

Here we give a brief explanation of the minimum-weight perfect matching (see Kelly et al. (2015) for a detailed explanation). We denote the measurement outcome of the $x$-th measurement qubit at the $y$-th cycle as $s_{x,y}$, and denote $m_{x,y}=s_{x,y}\oplus s_{x,y-1}$ with $s_{x,0}=0$, and $m_{\rm bnd}=\oplus m_{x,y}$. The minimum-weight perfect matching can then be regarded as finding the most probable error pattern in the following empirical model associated with a weighted graph $G$ on the set of vertices $V=\{(x,y)\}\cup\{{\rm bnd}\}$: starting with $m_{v}=0$ for all $v\in V$, an error occurs on each edge $(v,v^{\prime})$ independently with probability $e^{-w(v,v^{\prime})}$, where $w(v,v^{\prime})$ is the weight of $(v,v^{\prime})$. Whenever an error occurs on an edge $(v,v^{\prime})$, both $m_{v}$ and $m_{v^{\prime}}$ are flipped. For incoherent noises, the statistics of the actual circuit exactly follows the above model for an appropriate choice of $G$. In the case of the phenomenological model, $G$ is a uniformly weighted square lattice. For coherent noises, the statistics of the actual circuit does not exactly follow the simple model, and hence we need to choose $G$ heuristically. For simplicity, we choose the uniformly weighted square lattice for all the results shown in the main text. The possibility of using other graphs is discussed in Appendix C.

*Time efficiency

Since there are $O(n^{2})$ noise maps and syndrome measurements, and each of them takes at most $O(n^{3})$ steps, this scheme requires $O(n^{5})$ for each sampling. The time for simulation per sample with single thread of Intel Core i7 6700 takes about 20 $\rm ms$ with the parameter $(n,p,c)=(15,0.03,0)$ in the circuit-based model. For minimum-weight perfect matching, we used a library known as Kolmogorov’s implementation of Edmonds’ algorithm for minimum-weight perfect matching Edmonds (1965); Kolmogorov (2009).

Appendix B Appendix B: Obtaining the error threshold

For each model and each value of $c$, the error threshold $p_{\rm th}$ was determined by the following procedure. For various values of physical error probability $p$ and code size $n$, the logical error probability $p_{\rm L}(p,n)$ was computed as an average over 50,000 samples. Then it was fitted to the function $p_{\rm L}(p,n)=a+b(p-p_{\rm th})n^{1/d}$ Wang et al. (2003); Stephens (2014) around the error threshold, where $a,b,d$ and $p_{\rm th}$ are fitting parameters.

Appendix C Appendix C: Uniform and non-uniform weight in minimum-weight perfect matching

As mentioned in Process 3 of Appendix A, it is possible to find a weighted graph that faithfully reproduces the statistics of the actual circuit if the noises are incoherent. In the case of the circuit-based model, such a faithful graph has non-uniform weights and diagonal edges Kelly et al. (2015).

We constructed the faithful graph for the circuit-based model under incoherent noise, where each value of $e^{-w(v,v^{\prime})}$ is approximated with its leading $O(p)$ term for simplicity. We applied the decoder based on the constructed graph to the circuit-based model with various values of coherence $c$. The obtained error thresholds are shown in Fig. S1, together with the thresholds for the uniform-weight decoder. Compared to the error threshold using the uniform weight, the error threshold is improved for arbitrary values of $c$, but the amount of the improvement is small and the dependence of the error threshold on the noise coherence $c$ is also similar.

Appendix D Appendix D: Definition and calculation of the effective physical error probability

We define an effective physical error probability of a data qubit per cycle $p_{\rm eff}(p,c)$ as a marginal probability with which results of two measurement qubits neighboring a data qubit are flipped at a certain cycle from the results of the previous cycle. More precisely, we define $p_{\rm eff}(p,c)$ as the marginal probability of $m_{x,y}=m_{x+1,y}=1$, using the notation introduced in Process 3 of Appendix A. While $p_{\rm eff}(p,c)$ may depend on the values of $x$ and $y$, its leading term for small $p$ is independent of $x$ and $y$ (except $y=1$) and of the system size $n$. This leading term can be analytically obtained as $\frac{8}{3}(1+\frac{11}{6}c^{2})p$, and thus $\alpha=11/6$. This definition can be simply generalized to the case of the surface code, and we can analytically obtain the effective physical error probability of a data qubit per cycle since only a few noise maps and qubits are relevant to the leading term.

Appendix E Appendix E: Logical error scaling below threshold

If the noise is incoherent and the physical error probability $p$ is below the threshold, the logical error probability is expected to decay with code distance $d$ as

[TABLE]

where $\alpha$ is a polynomial function of $d$ Jones et al. (2012); Raussendorf et al. (2007). In this appendix, we numerically investigate whether this approximation is still valid for the coherent noises. Since $\alpha$ depends only weakly on $d$, the relation (S14) can be rewritten as

[TABLE]

We determined the above parameter $\lambda(d,p)$ from the computed values of $p_{\rm L}(d,p)$ and compared to the threshold $p_{\rm th}$ determined from the scaling ansatz. The result is shown in Fig. S2 for the case of the coherent noise ($c=1$) and the incoherent noise ($c=0$).

In both the phenomenological and the circuit-based model, Eq. (S14) is satisfied with the same level of approximation regardless of the degree of coherence in the noises.

Appendix F Appendix F: The non-unitary free-fermionic formalism of the surface code

Here we extend our scheme of the efficient classical simulation using non-unitary free-fermionic formalism to the surface code with a specific noise model including coherent errors. We consider the surface code with $d^{2}$ data qubits and $(d^{2}-1)$ measurement qubits for an odd number $d$. An example with $d=5$ is shown in Fig. S3(a).

The data qubits are located on the $d^{2}$ vertices of the colored faces. Each colored face corresponds to a measurement qubit. The measurement operator of each face is the product of $Y$ or $Z$ Pauli operators on the data qubits on its vertices. We denote the stabilizer operator of the blue face at the $r$-th row as $S_{r}^{\rm b}$ ($r=1,\ldots,d-1$). We denote the stabilizer operators of the red faces in the $c$-th column ($c=1,\ldots,d-1$) as $\{S^{{\rm r},c}_{i}\}$ ($i=1,\ldots,d-1$), and that of the green face as $S^{{\rm g},c}$. The rules for assigning the index $i$ will be explained later. For each cycle of syndrome measurement, the stabilizer operators as observables are measured via controlled gates and $Z$-basis measurement on the measurement qubits.

We consider an error model in which the code truly works as a fully quantum code, namely, including both $X$- and $Z$-type errors. More specifically, we assume the following phenomenological noise model. For each cycle, one of the following three types of errors occurs on each data qubit probabilistically, the Pauli $Y$ error, the Pauli $Z$ error, and an $X$-type coherent error as in Eq. (1) in the main text. We assume that the measurement qubits also suffer from the three types of errors probabilistically, the Pauli $X$, $Y$, and $Z$ errors. We also assume, for simplicity, that the syndrome measurement in the final cycle is error-free. This noise model can be considered as the phenomenological model in the main text with added Pauli $Y$ and $Z$ noises on both types of qubits, while limiting the $X$-type coherent errors to the data qubits. Apparently, such an error model including the coherent noise cannot be treated with the method based on the Gottesman-Knill theorem.

In contrast to the 1D repetition code in the main text where the logical error probability is the only parameter of interest, a circuit of fully quantum correction may be characterized in many different ways. Since we assumed that the final syndrome is correct, the whole circuit including a recovery operation can be viewed as a one-qubit channel $\mathcal{C}$ on the logical qubit space spanned by $\{\ket{0_{\rm L}},\ket{1_{\rm L}}\}$. The most general characterization of $\mathcal{C}$ is achieved if we learn the density matrix $\rho(\mathcal{C})$ on the logical qubit and an auxiliary qubit defined by

[TABLE]

with

[TABLE]

where ${\rm Id}$ is the identity channel. Any input-output relation can be calculated from $\rho(\mathcal{C})$, as well as the parameters such as various flip errors and fidelities.

Since the fermionic representation depends on the order of the data qubits, we assign numbers $1,2,\ldots,d^{2}(=n)$ to the data qubits in the order as shown in the gray boxes in Fig. S3(a). We also assume that the ancillary qubit is the $(n+1)$-th data qubit. Recall that the Majorana fermionic operators can be chosen as

[TABLE]

where $i=1,\ldots,n+1$. While the stabilizer operators for the blue and green faces are quadratic, those for the red faces are not. Our idea is to introduce a new equivalent set of stabilizer operators which are all quadratic. Let us consider the stabilizer operators in the $c$-th column, $\{S^{{\rm r},c}_{1},\ldots S^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$. With an appropriate rotation, each column can be considered as a part shown in Fig. S3(b). We assign the index of the stabilizer operators of the red faces in the column as shown in the figure. Let us introduce new red stabilizer operators defined from $\{S^{{\rm r},c}_{1},\ldots S^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$ as

[TABLE]

for $i=1,\ldots,d-1$. Since they are explicitly written in the form

[TABLE]

where $W$ is either $Z$ or $Y$, they are all quadratic in Majorana fermionic operators. The set $\{S^{{\rm r},c}_{1},\ldots S^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$ can be conversely expressed in terms of $\{\tilde{S}^{{\rm r},c}_{1},\ldots\tilde{S}^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$ as

[TABLE]

This indicates that $\{S^{{\rm r},c}_{1},\ldots S^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$ and $\{\tilde{S}^{{\rm r},c}_{1},\ldots\tilde{S}^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$ are equivalent sets as stabilizer generators. To be precise, let us introduce the projector onto the eigenspace of a stabilizer operator $S$ associated with syndrome bit $s$ by

[TABLE]

Then, it follows that

[TABLE]

if

[TABLE]

This implies that the syndrome measurement of $\{S^{{\rm r},c}_{1},\ldots S^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$ can be equivalently done by those of $\{\tilde{S}^{{\rm r},c}_{1},\ldots\tilde{S}^{{\rm r},c}_{d-1},S^{{\rm g},c}\}$ followed by calculation according to Eq. (S25).

Now we will show that the quantum error correction circuit can be efficiently simulated to compute $\rho(\mathcal{C})$. First, we show that the initial state $\ket{\phi_{\rm init}}$ is an FGS. The projection operator to the code space is given by

[TABLE]

which is an FGO. We consider the following logical $Z$ and $Y$ operators,

[TABLE]

and choose $\ket{0_{\rm L}}$ and $\ket{1_{\rm L}}$ to satisfy $L_{Z}\ket{0_{\rm L}}=\ket{0_{\rm L}}$, $L_{Z}\ket{1_{\rm L}}=-\ket{1_{\rm L}}$, and $L_{Y}\ket{0_{\rm L}}=i\ket{1_{\rm L}}$. Since $\ket{\phi_{\rm init}}$ is stabilized by $L_{Z}Z_{n+1}$, $-L_{Y}Y_{n+1}$, and the $d^{2}-1$ stabilizer operators, we have

[TABLE]

Since $L_{Z}Z_{n+1}$ and $L_{Y}Y_{n+1}$ are quadratic, $\ket{\phi_{\rm init}}\bra{\phi_{\rm init}}$ is an FGO, and hence $\ket{\phi_{\rm init}}$ is an FGS.

Simulation of the execution of the circuit is done as follows. The error on each data qubit is sampled, and if it is $X$-type, the state can be directly updated to a new FGS since the $X$-type error operator is an FGO. Pauli $Z$ and $Y$ operator is not an FGO, and hence it cannot be directly applied. When we implement $Z$ error on the $i$-th qubit, we apply $Z_{i}Z_{n+1}$ instead, which is an FGO, and the updated state can be calculated. We apply $Y_{i}Z_{n+1}$ for $Y$ error. The appearance of $Z_{n+1}$ should be recorded, and will be compensated in the final stage. A syndrome measurement of a stabilizer operator $S$ on a current FGS $\ket{\phi}$ without any error is done by calculating

[TABLE]

and sampling $s^{\prime}$ accordingly, and then calculating the updated state

[TABLE]

As we have seen, we are allowed to use $\tilde{S}^{{\rm r},c}_{i}$ instead of $S^{{\rm r},c}_{i}$, and to compute the syndrome for $\{S^{{\rm r},c}_{i}\}$ through Eq. (S25), which assures that $P(S,s^{\prime})$ is always an FGO. As for the Pauli error on the measurement qubit, it can be equivalently translated to a composition of the following operations: a bit-flip of the computed syndrome (for $X$ and $Y$), and Pauli errors on the data qubit (for $Y$ and $Z$). After all the cycles are executed, the recovery operator $R$ is calculated using the sampled syndrome values. $R$ can be written as a product of Pauli operators and hence computed by possible inclusion of $Z_{n+1}$.

After a single run of simulation, we obtain a final FGS $\ket{\phi_{\rm final}}$ and the record of the number $w$ of applications of $Z_{n+1}$. If $w$ is even, the FGS is an accurate sample of the desired state $\rho(\mathcal{C})$. If $w$ is odd, $Z_{n+1}\ket{\phi_{\rm final}}$ is an accurate sample of $\rho(\mathcal{C})$. Hence, all we need is to convert the description of $\ket{\phi_{\rm final}}$ as an $(n+1)$-qubit FGS to that as a two-qubit density operator. We choose the Pauli operators for the logical qubit as $L_{Z}$ and $L_{Y}$ defined in Eqs. (S27) and (S28), $L_{I}:=I^{\otimes n}$, and

[TABLE]

The density operator is then decomposed as

[TABLE]

where $W,W^{\prime}\in\{I,X,Y,Z\}$ and $A_{I,I}=1$. The coefficient $A_{X,X}$ can be obtained as follows. The final state $\ket{\phi_{\rm final}}$ is a +1 eigenstate of all the stabilizer operators. The operator $L_{X}$ can be transformed to $X^{\otimes n}$ by multiplying all the $(d-1)$ blue stabilizers, all the $(d-1)$ green stabilizers, and the $(d-1)^{2}/2$ red stabilizers composed of Pauli $Z$s. We thus have $L_{X}X_{n+1}\ket{\phi_{\rm final}}=X^{\otimes(n+1)}\ket{\phi_{\rm final}}$. Since $X^{\otimes(n+1)}$ commutes with all the FGOs and $\ket{\phi_{\rm init}}$ is a +1 eigenstate of $X^{\otimes(n+1)}$, we have $X^{\otimes(n+1)}\ket{\phi_{\rm final}}=\ket{\phi_{\rm final}}$, and hence $A_{X,X}=1$. This also implies $A_{W,W^{\prime}}=0$ if $L_{W}W^{\prime}_{n+1}$ anti-commutes with $L_{X}X_{n+1}$. The remaining 6 non-trivial coefficients are expectation values of $L_{X},X_{n+1},L_{Y}Y_{n+1},L_{Y}Z_{n+1},L_{Z}Y_{n+1}$, and $L_{Z},Z_{n+1}$ which are all expectation values of FGOs. These values can be calculated from the FGS description of $\ket{\phi_{\rm final}}$. We thus obtain the density operator $\ket{\phi_{\rm final}}\bra{\phi_{\rm final}}$.

An accurate sample of $\rho(\mathcal{C})$ is then given by applying the correction $Z_{n+1}^{w}$. The $\rho(\mathcal{C})$ is calculated as

[TABLE]

where $\rm{ave}(\cdot)$ represents averaging function over the samples.

Various parameters associated with channel $\mathcal{C}$ can be calculated from $\rho(\mathcal{C})$. For example, the entanglement fidelity $F(\mathcal{C})$ represents how the channel preserves the input quantum state. This is computed as

[TABLE]

Finally, we would like to introduce another fully quantum error correcting circuit which is efficiently simulated by our method. The circuit is a modification of the previous circuit. It uses the same surface code, but each of the new red syndromes $\{\tilde{S}^{{\rm r},c}_{i}\}$ is actually measured through a measurement qubit, instead of $\{S^{{\rm r},c}_{i}\}$. For this circuit, we can allow $X$-type coherent errors on measurement qubits, just as in the phenomenological model in the main text. In this case, the $X$-type coherent error can be absorbed by replacing the measurement operator $P(S,s)$ with

[TABLE]

where $\theta$ is a rotation angle dictated by the $X$-type coherent error. Since all the stabilizer operators are quadratic in the modified circuit, all of these operators are FGOs.

The second example shows that the 1D repetition code in the main text can be extended to a fully quantum code with no compromise on the noise model. Hence the applicability of our method relies neither on the simple structure of the repetition code nor on the absence of $Z$ and $Y$ errors. On the other hand, the comparison between the two examples of the surface code reveals an interesting trade-off. In the second example, the new syndrome measurements are non-local in the column direction, and only local in the row direction. In other words, it may be regarded as a 1D circuit with local stabilizer measurements. The first example is a true 2D circuit, but the efficient simulation seems to be possible only when the measurement qubits suffer no coherent errors. It is an open problem whether we can achieve both of them, efficient simulation of a truly 2D quantum error correction circuit with local syndrome measurement under coherent errors on both the data and measurement qubits.