Correcting coherent errors with surface codes

Surface codes are building blocks of quantum computing platforms based on 2D arrays of qubits responsible for detecting and correcting errors. The error suppression achieved by the surface code is usually estimated by simulating toy noise models describing random Pauli errors. However, Pauli noise models fail to capture coherent processes such as systematic unitary errors caused by imperfect control pulses. Here we report the first large-scale simulation of quantum error correction protocols based on the surface code in the presence of coherent noise. We observe that the standard Pauli approximation provides an accurate estimate of the error threshold but underestimates the logical error rate in the sub-threshold regime. We find that for large code size the logical-level noise is well approximated by random Pauli errors even though the physical-level noise is coherent. Our work demonstrates that coherent effects do not significantly change the error correcting threshold of surface codes. This gives more confidence in the viability of the fault-tolerance architecture pursued by several experimental groups. Coherent effects are shown not to play a significant role in error correction with quantum surface codes. To build a quantum computer, the quantum bit (qubit) has to be protected from external noise and steps have to be taken to detect and correct for errors. Surface codes are a type of quantum code that can correct for such errors. However, the models used to study such codes often fail to capture quantum coherent processes, which could play an important role. By performing large-scale simulations, Robert König from Technical University of Munich and an international team of collaborators show that coherent effects do not significantly impact the error correction in surface codes, giving confidence in the viability of this approach for developing fault-tolerance quantum computing architectures.


I. INTRODUCTION
Recent years have witnessed major progress towards the demonstration of quantum error correction and reliable logical qubits [1][2][3][4][5].Topological quantum codes such as the surface code [6,7] are among the most attractive candidates for an experimental realization, as they can be implemented on a two-dimensional grid of qubits with local parity check operators.
It is believed that such codes can tolerate a high level of noise [8][9][10] which is comparable to what can be achieved in the latest experiments [5].The general confidence in the noise-resilience of topological codes primarily rests on considerations of Pauli noise -a simplified noise model where errors are Pauli operators X, Y, Z drawn at random from some distribution.An example is the case where each qubit j experiences noise described by the channel with suitable probabilities x , y , z .This kind of noise can be fully described by the stabilizer formalism [11].
In pioneering work, Dennis et al. [8] exploited this algebraic structure to establish the first analytical threshold estimates, see also [12].The effect of Pauli noise also is efficiently simulable thanks to the Gottesman-Knill theorem, providing numerical evidence for high error thresholds of topological codes [13].The efficient simulability property has recently been extended beyond Pauli noise to random Cliffords and Pauli-type projectors [14].While such algebraically defined noise models are attractive from a theoretical viewpoint, they often do not correspond to noise encountered in real-world setups.They are -in a sense -not quantum enough: they model probablistic processes where errors act randomly on subsets of qubits.Rather than being of such a probabilistic (or incoherent) nature, noise in a realistic device will often be coherent, i.e., unitary, and can involve small rotations acting everywhere.A typical situation where this arises is if e.g., frequencies of oscillator qubits are mis-aligned: this results in systematic unitary over-or underrotations.On a single-qubit level, this means that (1) should be replaced by noise of the form with a suitable unitary operator U j ∈ SU (2).Since such errors generally cannot be described within the stabilizer formalism, understanding their effect on a given quantum fault-tolerant scheme is a challenging problem.
Prior theoretical work indicates that the difference between coherent and incoherent errors could be significant.In particular, it was observed [15][16][17][18][19] that coherent errors can lead to large differences between average-case and worst case fidelity measures suggesting that a critical reassessment of commonly used benchmarking measures is necessary.This observation motivates the question of how much coherence is present in the effective logical-level noise [20,21] experienced by encoded qubits.Depending on whether or not the logical noise is coherent one may choose different metrics for quantifying performance of a given fault-tolerant scheme.Significant progress has been made towards understanding the structure of the logical noise for concatenated codes [20][21][22].However, these studies are not directly applicable to large topological codes such as those considered here.
Brute-force simulations of coherent noise in small codes were presented in [23][24][25][26] for Steane codes and surface codes with up to 17 qubits.Simulating coherent errors by brute force clearly requires time (and memory) exponential in the number of qubits n.For the surface code, Darmawan and Poulin [27] proposed an algorithm with a runtime exponential in n 1/2 based on tensor networks, and simulated systems with up to 153 qubits.This algorithm can handle arbitrary noise (including e.g., amplitude damping).Unfortunately, its formidable complexity prevents accurate estimation of error thresholds, e.g., for the systematic rotations considered here.In [28], threshold estimates for the 1D repetition code were obtained.To our knowledge, there are no analogous threshold estimates for topological codes subject to coherent noise.
Our setup.Here we show that the effect of coherent errors in surface codes can be studied by means of polynomial-time algorithms.Specifically, we consider coherent errors in the context of two central tasks associated with error correction, namely (A) fault-tolerant storage of quantum information.
(B) fault-tolerant preparation of a logical basis state.
We shall consider a particular version of the surface code proposed in Refs.[29,30].A distance-d surface code has one logical qubit and n = d 2 physical qubits located at sites of a square lattice of size d × d with open boundary conditions.The code has local stabilizers X ⊗4 , X ⊗2 or Z ⊗4 , Z ⊗2 associated with faces of the lattice as shown in Fig. 1.The stabilizer located on a face f will be denoted B f .Logical Pauli operators X L and Z L acting on the encoded qubit can be chosen as X ⊗d and Z ⊗d applied to the left and the top boundary of the lattice respectively.The two-dimensional logical subspace is spanned by n- To specify the problem (A), consider a logical state ψ L initially encoded by the surface code and a coherent error applies some (unknown) unitary operator U j to each qubit j.To diagnose and correct the error without disturbing the encoded state we adopt the standard protocol based on the syndrome measurement.It works by measuring the eigenvalue (syndrome) s f = ±1 of each stabilizer B f on the corrupted state U |ψ L and then applying a Pauli-type correction operator C s depending on the measured syndrome s = {s f } f .The correction C s is computed by a classical decoding algorithm (for example, one may choose C s as a minimumweight Pauli error consistent with s).We note that the syndrome s is a random variable with some probability distribution p(s) since the error U maps the initial logical state to a coherent superposition of states with different syndromes.In this paper we only consider noiseless syndrome measurements.Accordingly, we assume that the correction C s always returns the system to the logical subspace resulting in some final logical state |φ s .For this problem we restrict to Z-rotation errors, that is, we assume that U j = exp(iη j Z) for some (unknown) angles η j .The restriction to Z-rotations is dictated by the limitations of our simulation algorithm.Thus we shall model a fault-tolerant storage by the following process: (i) prepare an initial logical state |ψ L (ii) apply a coherent error n j=1 exp (iη j Z) to |ψ L (iii) measure the eigenvalues of the stabilizers {B f } f , resulting in a syndrome s = {s f } f .
(iv) apply a Pauli correction C s returning the system to the logical subspace in some final state |φ s .
To assess how close the final state |φ s and the initial state |ψ L are, we seek a polynomial-time classical algorithm A which takes as input |ψ L and the rotation angles η 1 , . . ., η n , samples a syndrome s from the distribution p(s) specified by the measurement (iii), and outputs s as well as the associated final state |φ s (e.g.specified by its Bloch vector).By sampling sufficiently many syndromes, one can learn how frequently and in which ways error correction may fail in the presence of coherent noise.
To specify the problem (B), assume first that we have access to noise-free qubits and operations.In this case, the following standard protocol [10] prepares the encoded stabilizer state |+ L (the +1 eigenstate of X L ): (i) prepare the initial product state |+ ⊗n .
(ii) measure the eigenvalues of the stabilizers {B f } f , resulting in a syndrome s = {s f }.
(iii) apply a Pauli correction C s returning the system to the logical subspace in some final state |φ s .
Using the fact the the initial state is a +1 eigenvector of X L one can easily check that |φ s = |+ L for all s, see [10].How does this protocol fare in the presence of coherent noise?Let us consider a model where the initial product state cannot be prepared with perfect accuracy, but rather is obtained from |+ ⊗n by applying some unwanted unitary operators to every qubit.Thus (i) is replaced by , where ψ j are arbitrary single-qubit pure states.
In this case the final state |φ s may deviate from the target state |+ L resulting in a logical error.To assess the performance of this protocol, we seek a polynomial time algorithm B which, on input ψ 1 , . . ., ψ n ∈ C 2 , outputs a random syndrome s sampled from the distribution p(s) specified by the measurement (ii) together with the final logical state |φ s .
Our results.We construct algorithms A and B accomplishing the simulation tasks specified above.The runtime of these algorithms scales as O(n2 ), where we measure complexity in terms of the number of additions, multiplications, and divisions on complex numbers that are required1 .Using these algorithms, we perform the first numerical study of large topological codes subject to coherent noise, performing simulations for surface codes with up to n = 2401 physical qubits, see Table I for a timing analysis.This shows that efficient classical simulation of these fault-tolerance processes is possible, and allows us to extract key characteristics of these codes in the limit of large system size.
We apply algorithm A to study the effect of coherent noise on storage in the surface code.We show that the syndrome probability distribution p(s) is independent of the initial logical state ψ L whereas the final logical state has the form for some logical rotation angle θ s ∈ [0, π) depending on the syndrome s.We use the quantity as a measure of the logical error rate.We will see that P L is the average diamond-norm distance between the conditional logical channel ρ → e iθsZ ρe −iθsZ and the identity channel.For numerical simulations we consider translation-invariant coherent noise of the form (e iθZ ) ⊗n , where θ ∈ [0, π) is the only noise parameter.The Pauli correction C s was computed using the standard minimum weight matching decoder [8,31] with constant weights independent of θ.We are interested in the error threshold, that is, the maximum value θ 0 such that for any θ < θ 0 the logical error rate P L goes to zero in the limit n → ∞.We find the numerical estimate 0.08π ≤ θ 0 ≤ 0.1π.
Our numerical experiments confirm that, as expected, the quantity P L decays exponentially in the code distance for values θ < θ 0 below the threshold.Surprisingly, the threshold estimate Eq. ( 5) agrees very well with the so-called Pauli twirl approximation [32,33] where coherent noise of the form N (ρ) = e iθZ ρe −iθZ is replaced by dephasing noise D(ρ) = (1 − )ρ + ZρZ, with = sin 2 θ.
For the latter the threshold error rate is around 0 ≈ 0.11, see Ref. [8].Solving the equation 0 = sin 2 (θ 0 ) for θ 0 yields θ 0 ≈ 0.10π, in agreement with Eq. (5).At the same time, we observe that the Pauli twirl approximation significantly underestimates P L in the sub-threshold regime, confirming that coherence of noise may have a profound effect on a given fault-tolerant scheme, as was previously observed in [22,27].
Algorithm A allows us to investigate the probability distribution of logical rotation angle θ s defined in Eq. (3).We find that for large code sizes, this distribution concentrates around the two points {0, π/2} which correspond to the logical Pauli-type errors {I, Z L }.To get a deeper insight into this phenomenon, we introduce and numerically study associated measures of "incoherence".Our findings support the general conjecture that in the limit of large code distances, coherent physical noise gets converted into incoherent logical-level noise.
We apply Algorithm B to study the effect of coherent noise on logical state preparation in the case when the initial product state has the form (exp (iϕX) exp (iθZ)|+ ) ⊗n Here the angles θ, ϕ ∈ [0, π) specify the action of a coherent error on the ideal initial state |+ .The ideal protocol corresponds to θ = 0. We define the logical error rate P L as the average trace-norm distance between the final logical state φ s and the target state |+ L , see Section VI B for details.Our numerical results indicate that the error threshold can be described by a single function θ 0 (ϕ) such that the logical error rate P L goes to zero in the limit n → ∞ for any 0 ≤ θ < θ 0 (ϕ) and P L is lower bounded by a positive constant for θ > θ 0 (ϕ).We find the numerical estimate for all ϕ.This indicates that the threshold funciton θ 0 (ϕ) has a very mild (if any) dependence on ϕ.We investigate the behavior of P L in more detail for ϕ = 0 and obtain a more refined estimate 0.13π ≤ θ 0 (0) ≤ 0.14π.The quantity P L is observed to decay exponentially in the code distance in the sub-threshold regime.Outline.The remainder of this paper is structured as follows.Section II provides a high-level overview of our simulation algorithms.In Section III, we describe a representation of the surface code in terms of Majorana fermions.In Sections IV,V, we give the classical simulation algorithms A and B and analyze their complexity.In Section VI we discuss our numerical results.We conclude in Section VII.Appendix A contains a proof of a technical lemma.We provide some background on Majorana fermions and fermionic linear optics in Appendix B.

II. METHODS
Our main tool is a fermionic representation of the surface code proposed by Kitaev [34] and Wen [29].It works by encoding each qubit of the surface code into four Majorana fermions in a way that simplifies the structure of the surface code stabilizers.The Kitaev-Wen representation has previously been used by Terhal et al. [35] to design fermionic Hamiltonians with topologically ordered ground states.Here we show that this representation is also well-suited for the design of efficient simulation algorithms.The fermionic version of the surface code will be described in terms of Majorana operators c 1 , . . ., c 4n that obey the standard commutation rules c † p = c p , c 2 p = I, and c p c q = −c q c p for p = q.We will show that the error correction protocols considered in this paper can be decomposed into a sequence of O(n) elementary gates from a gate set known as a fermionic linear optics (FLO), see Refs.[36,37].It includes the following operations: 1. Initialize a pair of Majorana modes p, q in a basis state |0 satisfying ic p c q |0 = |0 .
3. Apply the projector Λ = (I + ic p c q )/2.Compute the norm of the resulting state.
It is well-known that quantum circuits composed of FLO gates can be efficiently simulated classically [36][37][38][39].The simulation runtime scales as O(n) for gates of type (1,2) and as O(n 2 ) for gates of type (3).For completeness, we describe the requisite simulation algorithms in Appendix B. By exploiting the geometrically local structure of the surface code we shall be able to reduce the number of modes such that at any given time step the simulator only needs to keep track of O(n 1/2 ) modes.Accordingly, each FLO gate can be simulated in time at most O(n).
Since the total number of gates is O(n), the total simulation time scales as O(n 2 ).

III. FROM QUBITS TO MAJORANA FERMIONS
A single Majorana mode p is described by a hermitian operator c p satisfying c 2 p = I.Operators c p associated with different modes anti-commute, see Appendix B for formal definitions.A system of four Majorana modes c 1 , c 2 , c 3 , c 4 can be used to encode a qubit using a stabilizer code with a single stabilizer and logical Pauli operators We shall refer to this encoding as a C4-code.
Consider a surface code with n qubits on a square d×d lattice, where n = d 2 .It can described by a planar graph G = (V, E, F ) with a set of n vertices V , a set of 2n − 2 edges E, and a set of n − 1 faces F .Qubits are located at vertices u ∈ V and stabilizers B f are located at faces f ∈ F of G. Consider a system of 4n Majorana modes c 1 , . . ., c 4n distributed over edges and vertices of G as shown on Fig. 2.There are exactly two paired modes located near the endpoints of every edge e ∈ E and four unpaired modes c 1 , c 2 , c 3 , c 4 located near the corners of the lattice as shown on Fig. 2. The paired modes are labeled as c 5 , c 6 , . . ., c 4n in an arbitrary order.Let us orient the edges of G as shown on Fig. 2. Suppose e ∈ E is an edge connecting some pair of modes c p , c q such that c p is the tail of e and c q is the head of e, see Fig. 2. Define the link operator Note that L e is hermitian and all link operators pairwise commute.Furthermore, each L e commutes with the unpaired modes c 1 , c 2 , c 3 , c 4 .By construction, a small neighborhood of each vertex u ∈ V contains a cluster of four modes, see Fig. 3.We shall denote this cluster Γ u .Define a vertex stabilizer where a particular product order is chosen for each vertex as shown on Fig. 3. Since |Γ u | = 4 for all u and the subsets Γ u are pairwise disjoint, vertex stabilizers are hermitian and pairwise commuting.We shall consider S 1 , . . ., S n as stabilizers for n independent copies of the C4-code defined in Eqs.(7,8) such that each qubit of the surface code is encoded into its own C4-code.Let X u , Z u , and Y u = iX u Z u be the logical Pauli operators for the qubit located at a vertex u, see Eq. (8).By definition, each of these logical operators has the form ic p c q for some pair of Majorana modes p, q ∈ Γ u , see Fig. 3.The logical operators X u , Z u are indicated by small arrows on Fig. 5.
Let P be a Pauli operator acting on the surface code qubits.We shall say that a Majorana operator P is a C4encoding of P if P can be obtained from P by replacing each single-qubit Pauli operator X u , Y u , Z u by its logical counterpart X u , Y u , Z u and, possibly, multiplying by the stabilizer S u .Given a single-qubit state ψ, one can define several encoded versions of ψ using the surface code, the C4-code, and the surface code concatenated with ncopies of the C4-code.We shall denote these encoded states ψ L , ψ, and ψ L respectively.These notations are summarized in Table II  The desired fermionic representation of the surface code is established in Lemmas 1,2,3 below.Consider a face f ∈ F and let ∂f ⊆ E be the boundary of f .Lemma 1.Let B f be the surface code stabilizer located on a face f .Then a C4-encoding of B f can be chosen as We illustrate Eq. ( 11) on Fig. 4. The lemma shows that measuring the surface code syndrome can be reduced (after the C4-encoding) to measuring eigenvalues of pairwise commuting link operators L e .We shall see that under certain circumstances such measurements can be efficiently simulated classically.
Proof.Consider a face f such that B f is a Z-stabilizer.Then B f = u∈f Z u .Consider a vertex u ∈ f and the C4-code located at u.The corresponding logical-Z operator Z u = ic p c q can be chosen such that both modes c p , c q are located on the boundary of f , see Fig. 3. Thus B f is proportional to the product of all modes located on the boundary of f .The same is true about the operator By construction, the boundary ∂f alternates between link operators L e and logical-Z operators Z u , see Fig. 4 for an example.For each link operator L e with e ∈ ∂f define a quantity ω f (L e ) = ±1 such that ω f (L e ) = −1 iff e is oriented clockwise with respect to f .Likewise, for each logical-Z operator Z u = ic p c q lying on the boundary of f define a quantity ω f (Z u ) = ±1 such that ω f (Z u ) = −1 iff an arrow c p → c q is oriented clockwise with respect to f .See Fig. 5 for examples of such arrows.A simple computation shows that e∈∂f L e = −ω f B f , where Thus we need to check that ω f = −1 for each face f ∈ F .In other words, the boundary of each face must have an odd number of arrows oriented clockwise.Direct inspection shows that this is indeed the case for the distance-3 code, see Fig. 5.By translation invariance, this also holds for all code distances.The same arguments apply to Xtype stabilizers.
We shall need an analogue of Lemma 1 for logical operators of the surface code.Lemma 2. Let X L and Z L be the logical operators of the surface code located on the left and the top boundary.
Then C4-encodings of X L and Z L can be chosen as where LEFT and TOP are the subsets of edges lying on the left and the top boundaries of the lattice, see Fig. 6.

FIG. 6. The sets of edges LEFT (red) and TOP (blue).
Proof.Let us add a "logical edge" connecting the modes c 2 and c 3 to the graph G, see Fig. 2.This creates an extra "logical face" f attached to the top boundary of the lattice.The new edge carries a link operator L e = ic 2 c 3 .The same arguments as in the proof of Lemma 1 show that Z L = −ω f e∈∂f L e , where ω f = ±1 is the parity of the number of arrows lying on the boundary of f and oriented clockwise with respect to f .From Fig. 5 one gets ω f = −1.The same argument applies to X L .Suppose ψ is a single-qubit state.The following lemma allows one to switch between different encodings of ψ defined in Table II by measuring syndromes of the vertex stabilizers.Informally, it asserts that a qubit initially encoded into the four unpaired Majorana modes c 1 , c 2 , c 3 , c 4 at the corners of the lattice can be "injected" into the logical subspace of the surface code by tensoring in unentangled pairs of modes located on edges and measuring syndromes of the vertex stabilizers.Lemma 3. Let φ link be the state of 4n − 4 Majorana modes c 5 , c 6 , . . ., c 4n stabilized by all link operators, Let ψ be any single-qubit state.Then Here we used the notations from Table II.
The lemma will allow us to replace the initial logical state ψ L in the error correction protocol by a simpler state φ link at the cost of measuring certain additional stabilizers.We shall see that the state φ link is a fermionic Gaussian state, see Section B for details.Furthermore, a state obtained from φ link by applying a coherent error (encoded by the C4-code) is also Gaussian.These features will be instrumental for our simulation algorithm.
Proof.Suppose first that |ψ = |0 .Let us add a pair of "logical edges" connecting modes c 2 , c 3 and c 1 , c 4 to the graph G.This creates an extra pair of "logical faces" attached to the top and the bottom boundaries.The new edges carry link operators ic 2 c 3 and ic 1 c 4 .Lemmas 1,2 imply that the state |ψ ⊗|φ link is stabilized by operators B f and Z L .Furthermore, since these operators commute with all vertex stabilizers, the state on the right-hand side of Eq. ( 14) is stabilized by B f and Z L .Since it is also stabilized by all S u , this state has the same set of stabilizers as |0 L , which proves Eq. ( 14) for |ψ = |0 .Note that see Lemma 2. Furthermore, X L commutes with all vertex stabilizers.Thus applying X L to both sides of Eq. ( 14) with |ψ = |0 proves Eq. ( 14) for |ψ = |1 .By linearity, it holds for all ψ.

IV. LOGICAL STATE PREPARATION
We shall first describe the algorithm for simulating the logical state preparation because it is much simpler than the storage simulation.For each syndrome s = {s f } f ∈F define a syndrome projector It projects onto the subspace spanned by n-qubit states with the syndrome s.Note that s Π s = I.Let Π 0 be the projector onto the logical subspace of the surface code.Since the Pauli correction C s maps any state with a syndrome s to the logical subspace, it must satisfy Here and below we assume that C † s = C s .Suppose that our initial state has the product form Here ψ j are arbitrary single-qubit states.Our goal is to sample a syndrome s from the probability distribution and compute the final logical state conditioned on the syndrome.The latter has the form Let us first discuss how to simulate the syndrome measurement.We shall encode each qubit into the C4-code as discussed in Section III.Let |ψ a be the encoded version of |ψ a .Using the Euler angle decomposition one can write |ψ a = e −iαX e −iβZ e −iγX |0 .Replacing Pauli operators by their C4-encodings defined in Eq. (8) gives The state |0 is stabilized by Z = ic 2 c 3 and ZS = ic 1 c 4 , see Eq. (8).Thus the states |ψ a can be prepared using only FLO gates.Applying this independently to each qubit gives a sequence of O(n) FLO gates that prepares the state Clearly, measuring syndromes of stabilizers B f on the state |ψ is equivalent to measuring syndromes of the encoded stabilizers B f on the state |ψ .By Lemma 1, the latter can be reduced to measuring syndromes m e = ±1 of the link operators L e and then classically computing face syndromes s f = e∈∂f m e .Since each link operator has a form L e = ic p c q , measuring the eigenvalue of L e is a FLO gate.Thus we have realized the full syndrome measurement using O(n) FLO gates.
It remains to compute the final logical state φ s .Let b s = (b x s , b y s , b z s ) be the Bloch vector of φ s such that Below we focus on computing b x s .Define Here in the first equality we noted that X L commutes with Π s .In the second equality we encoded every qubit using the C4-code.Let m ∈ {+1, −1} E be the combined syndrome of link operators measured in the first part of the algorithm such that m e is the measured eigenvalue of L e .Define the corresponding syndrome projector We claim that Here LEFT is the subset of edges lying on the left boundary of the lattice, see Fig. 6.The ratio in Eq. ( 20) can be computed by taking the normalized state Π link m |ψ obtained after measuring the link syndromes (the first part of the algorithm) and measuring the eigenvalue of ic 1 c 2 .The latter requires a single FLO gate.This gives the desired value of b x s .Likewise, measuring the eigenvalues of ic 2 c 3 and −ic 1 c 3 on the final state Π link m |ψ gives the remaining components of the Bloch vector b z s and b y s respectively.
To conclude, FLO gates enable simulation of the syndrome measurement and computation of the final logical Bloch vector conditioned on the measured syndrome.
The resulting FLO circuit can be simulated classically in time O(n 3 ) since it includes O(n) projection gates (1/2)(I + m e L e ).The simulation runtime can be significantly improved using the following simple observations.First, once a link operator L e = ic p c q has been measured, the modes c p , c q are completely disentangled from the rest of the system.Such disentangled modes can be removed from the simulator reducing the total number of modes and the computational cost of subsequent steps.Second, we can exploit the fact that the initial state ψ has a product form.In particular, a four-mode cluster Γ u that supports the state ψ u only needs to be loaded into the simulator at a time step when some mode p ∈ Γ u participates in the measurement of a link operator.Thus at any given time step the simulator only needs to keep track of "active" clusters Γ u such that at least one mode p ∈ Γ u has been measured and at least one mode q ∈ Γ u has not been measured.One can easily choose the order of measurements such that the number of active clusters is O(n 1/2 ) at any time step (for example, one can measure link operators column by column).Accordingly, the cost of simulating a single FLO gate is at most O(n).Since the number of gates is O(n), the total simulation cost is O(n 2 ).
It remains to prove Eq. ( 20).Let δ be a map from link syndromes to the corresponding face syndromes, that is, Below we prove the following Proposition 1. Suppose m and m are link syndromes such that δ(m) = δ(m ).Then there exists a subset of vertices W ⊆ V such that We postpone the proof until the end of the section.Let m and m be any link syndromes with δ(m) = δ(m ).Then the proposition implies that since S u ψ = ψ for all u.Likewise, since X L commutes with T .From Eqs. (19,21,22,23) one gets for any link syndrome m such that δ(m) = s.In particular, one can choose m as the link syndrome measured in the first part of the algorithm.Substituting X L from Lemma 2 gives the desired result Eq. (20).It remains to prove Proposition 1.
Proof.By linearity, we can assume that m is the trivial syndrome, that is, m e = 1 for all e ∈ E. Then δ(m ) = δ(m) iff m is a flat connection on the surface code lattice, that is, m is a 1-chain with Z 2 coefficients such that the Z 2 -valued magnetic flux through every face is +1.Since the lattice is topologically trivial, any flat connection is a co-boundary of some 0-chain W ⊆ V .Since a vertex stabilizer S u flips the link syndromes on all edges incident to u, the condition that m is a co-boundary of W is equivalent to Π link m = T Π link m T for T = u∈W S u .

V. STORAGE OF A LOGICAL STATE
Assume now that our initial state has the form where ψ L is some (unknown) logical state of the surface code and η 1 , . . ., η n are arbitrary rotation angles.The unitary U describes a coherent error applied to each qubit before the syndrome measurement.Our goal is to sample a syndrome s from the probability distribution and compute the final logical state conditioned on the syndrome, Clearly, since ψ contains only Z-type errors, the observed syndrome of Z-stabilizers is always trivial.Accordingly, we shall assume that the correction C s is a Z-type Pauli.We shall need the following fact.
Lemma 4. The probability p(s) does not depend on the initial logical state ψ L .The map ψ L → φ s is a logical Z-rotation by some angle θ s ∈ [0, π), that is, We shall refer to θ s as a logical rotation angle.Here and below all states are defined modulo an overall phase factor.We defer the proof of the lemma to Appendix A.

A. Simulating the syndrome measurement
Let us first discuss how to sample s from p(s).By Lemma 4, p(s) does not depend on ψ L , so below we set |ψ L = |+ L .Since only syndromes of X-stabilizers may be non-trivial, it suffices to measure eigenvalues m u = ±1 of single-qubit Pauli operators X u and then classically compute the face syndrome s f = u∈f m u for each face f that supports an X-stabilizer.Let m ∈ {+1, −1} V be the combined X-measurement outcome and p x (m) be the probability of an outcome m.We have Let us order qubits column by column as shown on Fig. describing the first t qubits and let be the conditional distribution of m t given m 1 , . . ., m t−1 .
Let us show how to sample m t from the distribution Eq. ( 29) using FLO gates.Partition the set of all qubits as [n] = AB where A = {1, . . ., t} and B = [n] \ A. We shall write Then Encoding each qubit into the C4-code as discussed in Section III gives The Majorana representation of the logical state ψ L defined in Lemma 3 gives where γ is a normalizing coefficient depending only on n, φ link is a product state defined in Eq. ( 13), and ψ L is the basis state of modes c 1 , c 2 , c 3 , c 4 stabilized by ic 1 c 2 and ic 3 c 4 (recall that we have chosen |ψ L = |+ L ).The state |ψ L ⊗ |φ link can be prepared using FLO gates since it is stabilized by two-mode operators L e , ic 1 c 2 and ic 3 c 4 .
To simplify the notation, we shall absorb the pairs ic 1 c 2 and ic 3 c 4 into φ link .Accordingly, below we assume that φ link is a state of 4n modes defined as (32) Plugging Eq. ( 31) into Eq.(30) gives where Ω = u∈V 1 2 (I + S u ) is the projector onto the codespace of the C4-code.Write Ω = Ω A Ω B , where Since U A and Π x A commute with Ω A , one gets Assume first that B = ∅.Expand the projector Ω B as Here we used the fact that B is a connected set for any t, see Fig. 7.
Substituting the expansion of Ω B into Eq.( 34) and using the above observation gives where γ = 2 −|B| γ.Next we note that φ link is stabilized by S A S B since the latter coincides with the product of all link operators L e (including ic 1 c 2 and ic 3 c 4 ) and φ link is stabilized by any link operator.Thus one can replace S B by S A in Eq. (35).However, since S A commutes with U A and S A Ω A = Ω A , we arrive at Using the identity one finally gets where G 2a−1 , G 2a are operators acting non-trivially only on the subset of modes Γ a , namely and Noting that Z a , X a and X a S a have the form ic p c q for some p, q ∈ Γ a , see Eqs. (7,8), one concludes that Eq. ( 37) includes only FLO gates.
The above arguments also apply to the case B = ∅, that is, t = n.The only difference is that now Eq. ( 35) has no term S B and thus Eq. ( 37) becomes Let us discuss the cost of sampling m from p x (m).First, observe that any single-qubit marginal state of ψ L is maximally mixed since the surface code has no stabilizers of weight one.Since the same is true about the state U |ψ L , one can pick the first measurement outcome m 1 at random from the uniform distribution.Suppose we have already sampled m 1 , . . ., m t−1 and the simulator's current state is Plugging Eqs.(37,38) into Eq.( 29) gives where κ t = 2 for t < n and κ n = 1.The conditional probability of the outcome m t = 1 can be computed by simulating two more FLO gates G 2t , G 2t−1 starting from the state φ t−1 which takes time O(n 2 ).Once the conditional probability is computed, one can sample m t by tossing a suitably biased coin.This produces the next syndrome m t together with the next state φ t .After n iterations one gets the desired sample m from p x (m).Since the above algorithm uses O(n) FLO gates, the simulation runtime scales as O(n 3 ).As before, we can reduce the runtime to O(n 2 ) by exploiting the fact that the initial state φ link has a product form and by observing that once a pair of modes have been measured, it can be removed from the simulator.Indeed, consider some intermediate step t and let j be the column of the lattice that contains t-th qubit.Then all modes in the columns j + 2, . . ., d can be grouped into unentangled pairs located on edges such that each pair is stabilized by the respective link operator L e .Such unentangled pairs do not need to be loaded into the simulator.Likewise, all modes in the columns 1, . . ., j − 2 have already been measured and can be removed from the simulator.Thus at any given time step the number of "active" modes that needs to be simulated is only O(n 1/2 ).Accordingly, the cost of simulating a single FLO gate is at most O(n).Since the number of gates is O(n), the total simulation cost is O(n 2 ).
Remark 1: The same reasoning shows that the probability p x (m) of any given outcome m can be computed up to the normalizing coefficient γ in time O(n 2 ) by simulating the FLO circuit defined in Eq. (38).Furthermore, the normalizing coefficient γ depends only on n (but not on the rotation angles η a ).
Remark 2: Below we shall also use a slightly modified version of the above algorithm where the initial state ψ L is an eigenvector of the logical-Y operator.The modified version is exactly the same as above except that the initial state φ link in Eq. ( 32) is defined as (41) This corresponds to initializing the unpaired modes in the logical-Y state.

B. Computing the logical rotation angle
It remains to show how to compute the logical rotation angle θ s for a given syndrome s.Let us first initialize the logical qubit in the X-basis, that is, we choose |ψ L = |+ L .Let φ s be the final logical state defined in Eq. ( 26).Define logical amplitudes and Here we used Lemma 4. Then Using Eq. ( 26) and the identity Let U + = C s U and U − = Z L C s U .By definition, U ± are products of single-qubit Z rotations: Let us expand the logical state |+ L in the Z-basis: where L is the set of basis states x ∈ {0, 1} n that obey Zstabilizers of the surface code (we do not need an explicit formula for L).Since U ± is diagonal in the Z-basis, Substituting this into Eq.( 44) gives Since U ± is a tensor product of Z-rotations, p ± is a special case of the probability p x (m) defined in Eq. ( 28) with m u = 1 for all u.We have already shown that one can compute γ −1 p ± in time O(n The same arguments as above show that Since U ± is a product of Z-rotations, q ± is a special case of the probability p x (m) defined in Eq. ( 28) with m u = 1 for all u and the initial state ψ L chosen as a logical Ybasis state.We have already shown that one can compute q ± in time O(n 2 ) see Remarks 1,2 at the end of Section V A. Thus tan 2 (θ s − π/4) = q − /q + is computable in time O(n 2 ).Combining Eqs.(47,50) gives This determines the logical rotation angle modulo π.

VI. NUMERICAL RESULTS
We implemented the algorithms described above for translation-invariant coherent noise and surface codes with distance 5 ≤ d ≤ 49.The smallest distance d = 3 was skipped because of strong finite-size effects (note that the considered surface codes are only defined for odd values of d).We used the maximum distance d = 37 for storage simulations and d = 49 for state preparation simulations.The logical error rate P L was estimated by the Monte Carlo method with at least 50, 000 syndrome samples per data point.(The only exception is Fig. 12 where we used 5, 000 syndrome samples per data point.)

A. Numerical results for storage
Consider first the protocol A for storage of a logical state in the presence of Z-rotation errors.In this section we only consider translation-invariant errors of the form exp (iθZ) ⊗n , where θ is the only noise parameter.Recall that we define the logical error rate as where p(s) is the probability of observing a syndrome s and θ s is the logical rotation angle conditioned on the 0.06 0.08 0.10 0.12 0.14 Z-rotation angle / that describes the residual logical error for a given syndrome s.Let • denote the diamond-norm [40] on the space of quantum channels and id be the single-qubit identity channel.The identity Λ s − id = 2| sin θ s | shows that P L coincides with the average diamond-norm distance between the conditional logical channel and the identity channel, Using the symmetries of the surface code one can easily check that P L is invariant under flipping the sign of θ.Accordingly, it suffices to simulate θ ≥ 0.
Our numerical results for the logical error rate are presented in Fig. 8.The data suggests that the quantity P L decays exponentially in the code distance d for θ < θ 0 , where 0.08π ≤ θ 0 ≤ 0.1π (54) can be viewed as an error correction threshold.We observe the exponential decay of P L as a function of d in the sub-threshold regime.Although the logical error rate P L is a meaningful measure of how well the initial logical state is preserved, it provides no insight into the structure of residual logical errors.Algorithm A gives us a unique opportunity to investigate the logical-level noise since it outputs both the syndrome s and the the logical rotation angle θ s conditioned on s.Fig. 9 shows the empirical probability distribution of θ s obtained by sampling 10 6 syndromes s for the physical Z-rotation angle θ = 0.08π (which we expect to be slightly below the threshold).We compare the cases d = 9 and 25.In both cases the distribution has a sharp peak at θ s = 0 (equivalent to θ s = π).This peak indicates that error correction almost always succeeds in the considered regime.For ease of visualization, we truncated the peak at θ s = 0 on the histograms.It can be seen that increasing the code distance has a dramatic effect on the distribution of θ s .The distance-9 code has a broad distribution of θ s meaning that the logical-level noise retains a strong coherence.On the other hand, the distance-25 code has a sharply peaked distribution of θ s with a peak at θ s = π/2 which corresponds to the logical Pauli error Z L .Such errors are likely to be caused by "ambiguous" syndromes s for which the minimum weight matching decoder makes a wrong choice of the Pauli correction C s .We conclude that as the code distance increases, the logical-level noise can be well approximated by random Pauli errors even though the physical-level noise is coherent.
To investigate this effect more systematically, it is desirable to have a metric quantifying the degree of coherence present in the logical-level noise.To this end let us consider the twirled version of the logical channel Λ s , and the corresponding logical error rate FIG.10.Coherence ratio P L /P L twirl for the conditional logical channel (left) and for the average logical channel (right).In both cases increasing the code distance makes the logical-level noise less coherent.on {0, π/2}, that is, when the logical noise is incoherent.It is therefore natural to measure coherence of the logical noise by the ratio P L /P L twirl .This "coherence ratio" is plotted as a function of θ on Fig. 10(a).The data indicates that the coherence ratio decreases for increasing system size approaching one for large code distances.This further supports the conclusion that the logical noise has a negligible coherence.Finally, in Fig. 10(b), we show the analogous quantity for the average logical noise chan-nel [20] defined as Λ = s p(s)Λ s .
This average channel provides an appropriate model for the logical-level noise if the environment has no access to the measured syndrome.This may be relevant, for instance, in the quantum communication settings where noise acts only during transmission of information.Thus In both cases we used the minimum weight matching decoder.The plot demonstrates that applying the Pauli twirl approximation to the physical noise significantly underestimates the logical error rate in the sub-threshold regime.
one can alternatively define the coherence ratio as where Λ twirl is the Pauli-twirled version of Λ.In our case Λ(ρ) = (1 − )ρ + ZρZ + iδ(Zρ − ρZ), where = s p(s) sin 2 (θ s ) and δ = s p(s) sin (2θ s )/2, see Eq. (53).A simple calculation yields The coherence ratio of the average logical channel is plotted as a function θ on Fig. 10(b).It provides a particularly strong evidence that in the limit of large code distances, coherent physical noise gets converted into incoherent logical noise.
Finally, let us compare logical error rates P L computed for coherent physical noise N (ρ) = e iθZ ρe −iθZ and its Pauli-twirled version N twirl (ρ) = (1 − )ρ + ZρZ with = sin 2 (θ).Applying the Pauli twirl at the physical level amounts to ignoring the coherent part of the noise.Let P L (N twirl ) be the logical error rate corresponding to N twirl .The plot of P L (N twirl ) and the ratio P L /P L (N twirl ) are shown on Fig. 11.It can be seen that applying the Pauli twirl approximation to the physical noise gives an accurate estimate of the error threshold but significantly underestimates the logical error probability in the subthreshold regime.We conclude that coherence of noise may have a profound effect on the performance of large surface codes in the sub-threshold regime which is particularly important for quantum fault-tolerance.

B. Numerical results for state preparation
Next consider the protocol B for preparing the logical basis state |+ L by performing syndrome measurements on the initial product state (exp (iϕX) exp (iθZ)|+ ) ⊗n .
Here ϕ, θ ∈ [0, π) are noise parameters.The ideal protocol corresponds to θ = 0. Define the logical error rate P L as the average trace-norm distance between the final logical state |φ s φ s | and the target state |+ L + L |. Equivalently, Since the considered noise model generates correlations between X-and Z-syndromes, we opted not to use the minimum-weight matching decoder (which treats X-and Z-syndromes independently).Instead, we used a simplified decoder that chooses a Pauli correction C s such that the final logical state φ s always obeys φ s |X L |φ s ≥ 0. The simplified decoder is optimal in the sense that it minimizes the logical error rate P L under the constraint that C s is a Pauli operator.Note that 0 ≤ P L ≤ √ 2 for all noise parameters.The symmetries of the surface code imply that P L , considered as a function of θ and ϕ, is invariant under transformations Thus it suffices to simulate the region 0 ≤ θ, ϕ ≤ π/4.Our numerical results for a fixed code distance d = 39 are presented on Fig. 12.The data supports a natural conjecture that the asymptotic behavior of P L in the limit n → ∞ can be characterized by a single threshold function θ 0 (ϕ) such that lim n→∞ P L = 0 for θ < θ 0 (ϕ) and P L is lower bounded by a positive constant independent of n for θ > θ 0 (ϕ), see also Fig. 13.From Fig. 12 one gets an estimate 0.1π ≤ θ 0 (ϕ) ≤ 0.15π for all ϕ.The data indicates that θ 0 (ϕ) has a very mild (if any) dependence on ϕ.To test the above conjecture, we performed more detailed simulations for ϕ = 0, see Figs. 13,14, obtaining a more refined estimate 0.13π ≤ θ 0 (0) ≤ 0.14π.

VII. CONCLUSIONS
Our work extends the range of noise models efficiently simulable on a classical computer.It allows -for the first time -to numerically investigate the effect of coherent errors in the regime of large code sizes which is important for reliable error threshold estimates.Our simulation algorithms make no assumptions about the particular decoder used.Hence the proposed approach should be universally applicable to benchmarking the performance of different fault-tolerance strategies in the presence of coherent noise.
Our numerical results spell good news for quantum engineers pursuing surface code realizations: thresholds for state preparation and storage are reasonably high, suggesting that coherent noise is not as detrimental as one could expect from the previous studies.The numerical investigation of the logical-level noise gives rise to a conceptually appealing conjecture: error correction converts coherent physical noise to incoherent logical noise (for large code sizes).Whether this is an artifact of the considered error correction scheme or manifestation of a more general phenomenon is an interesting open question.
Although we simulated only translation-invariant noise models, all our algorithms apply to more general qubitdependent noise.This enables numerical study of recently proposed state injection protocols [41], e.g.preparation of logical magic states, in the presence of coherent errors.Another possible application could be testing the so-called disorder assisted error correction method [42][43][44] where artificial randomness introduced in the code parameters suppresses coherent propagation of errors due to the Anderson localization phenomenon.We leave as an open question whether our algorithms can be extended to more general coherent noise models such as those including systematic cross-talk errors.The above facts are sufficient to simulate any quantum circuit composed of FLO gates, as defined in Section II, and provide all details necessary for implementation of our algorithms A and B.

FIG. 1 .
FIG. 1. Surface codes with distance d = 3 and d = 5.Qubits and stabilizers are located at sites and faces respectively.A stabilizer B f located on a face f applies X (black faces) or Z (white faces) to each qubit on the boundary of f .Logical Pauli operators XL (red) and ZL (blue) have support on the left and the top boundary.

FIG. 2 .
FIG.2.Solid circles represent paired (blue) and unpaired (red) Majorana modes.An edge oriented from cp to cq defines a link operator Le = icpcq.The pattern extends to larger codes in a translation invariant fashion.

FIG. 4 .
FIG.4.Left: C4-encoding of a Z-type surface code stabilizer B f .Right: The same operator represented as a product of link operators Le over the boundary of f .

FIG. 8 .
FIG.8.Logical error rate P L for storage of quantum states.We consider distance-d surface codes subject to coherent errors exp (iθZ) on each qubit.

FIG. 9 .
FIG. 9.These histograms show the empirical probability distribution of logical rotation angles θs for the code distance d = 9 (left) and d = 25 (right).The histograms use the same noise parameter θ = 0.08π.For ease of visualization, we truncated the main peak at θs = 0.

FIG. 11 .
FIG. 11.Comparison between the logical error rates P L and P L (N twirl ) computed for coherent noise N (ρ) = e iθZ ρe −iθZ and its Pauli twirled version N twirl (ρ) = (1 − )ρ + ZρZ with = sin 2 (θ).In both cases we used the minimum weight matching decoder.The plot demonstrates that applying the Pauli twirl approximation to the physical noise significantly underestimates the logical error rate in the sub-threshold regime.

4 FIG. 12 .
FIG.12.Preparation of the logical basis state |+L by performing syndrome measurements on the initial product state (exp (iϕX) exp (iθZ)|+ ) ⊗n .The color represents the logical error rate P L defined in Eq. (57).The dark blue and red regions represent a "fault-tolerant" phase where the final logical state is close to |+L and |0L respectively.Here we consider a fixed code distance d = 39.The plot was generated by sampling 5, 000 syndromes for each pair (θ, ϕ).

TABLE II
. Encoded versions of a single-qubit state ψ.
Consider a link operator L e associated with some edge e such that both endpoints of e belong to B. Such a link operator commutes with all operators acting on A. Note that L e commutes with all vertex stabilizers S u except for the two stabilizers located at the endpoints of e.It follows that L e S C L e = −S C if exactly one end-point of e belongs to C. Furthermore, since φ link is stabilized by all link operators one infers that φ link |O A S C |φ link = 0 for any operator u∈CS u .