(y Selected for a Viewpoint in Physics PHYSICAL REVIEW X 2, 021004 (2012)

Strong Resilience of Topological Codes to Depolarization

H. Bombin,1 Ruben S. Andrist,2 Masayuki Ohzeki,34 Helmut G. Katzgraber,52 and M. A. Martin-Delgado6

1 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada 2Theoretische Physik, ETH Zurich, CH-8093 Zurich, Switzerland 3Department of Systems Science, Graduate School of Informatics, Kyoto University, Yoshida-Honmachi,

Sakyo-ku, Kyoto 606-8501, Japan 4Dipartimento di Fisica, Universitci di Roma ''La Sapienza'', P.le Aldo Moro 2, 00185, Roma, Italy 5Department of Physics and Astronomy, Texas A&M University, College Station, Texas 77843-4242, USA 6Departamento de Física Teórica I, Universidad Complutense, 28040 Madrid, Spain (Received 8 February 2012; published 30 April 2012)

The inevitable presence of decoherence effects in systems suitable for quantum computation necessitates effective error-correction schemes to protect information from noise. We compute the stability of the toric code to depolarization by mapping the quantum problem onto a classical disordered eight-vertex Ising model. By studying the stability of the related ferromagnetic phase via both large-scale Monte Carlo simulations and the duality method, we are able to demonstrate an increased error threshold of 18.9(3)% when noise correlations are taken into account. Remarkably, this result agrees within error bars with the result for a different class of codes—topological color codes—where the mapping yields interesting new types of interacting eight-vertex models.

DOI: 10.1103/PhysRevX.2.021004Subject Areas: Computational Physics, Quantum Information, Statistical Physics

I. INTRODUCTION

Moore's law has accurately described the speedup of current computer technologies for half a century, yet this speedup is slowly coming to an end due to transistor limitations. A promising alternative is given by quantum computers. However, the qubit manipulations required for information processing and communication are prone to errors because qubits are more sensitive to noise than their classical counterparts. Consequently, protecting qubits has become an issue of paramount importance for the success of quantum computation. The effects of single-qubit operations can be decomposed into three processes—bit flips, phase flips, as well as a combination thereof, which can be represented by the three Pauli matrices ax, az, and ay, respectively. This is in contrast to classical bits, which can suffer only from a single type of error, namely, bit flips.

More generally, the notion of a noisy channel is instrumental in characterizing the disturbing effects on physical qubits. Such a quantum channel can be described by specifying the probability (or "qubit error rate'') p for each of the aforementioned noise types. For instance, if only ax occurs, then we have a bit-flip channel. In this paper, we are interested in channels of the form

D p (p) = (1 - p)p + X

pwawpaw,

Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

where the density matrix p fully describes the quantum state, and the probability for each type of error to occur is pw 2 [0, 1] with p := px + py + pz. The depolarizing channel exhibits equal probability pw = p/3 for each error type to arise. That being so, the depolarizing channel is more general than the bit-flip channel, because it allows for the unified, correlated effect of all three basic types of errors. The implications of this error model for the performance of a quantum code remains an open problem. In addition, the depolarizing channel plays a fundamental role in quantum-information protocols where noise has to be taken into account, including quantum cryptography [1,2], quantum distillation of entanglement [3], and even quantum teleportation [4]. Experimentally, controllable depolarization has been realized recently in photonic quantum-information channels [5] with a Ti:sapphire pulsed laser and nonlinear crystals, as well as 2-qubit Bell states [6]. Here we compute the effects of depolarization on a set of entangled qubits.

A. Topological codes

The goal of quantum error correction [7,8] is to protect quantum information from decoherence. One approach using topology is based on encoding (few) logical qubits in a particular state subspace of (typically many) physical qubits which is not disturbed directly by noise. Such a suitable subspace of states can be defined in terms of a set of commuting observables, called check operators,

S¿ = aV2

each being a projective measurement with respect to the code subspace (i.e., the eigenvalue signals errors on participating qubits). Investigating all stabilizers S' allows one

w=x,y,z

2160-3308/12/2(2)/021004(10)

021004-1

Published by the American Physical Society

to limit the set of possible errors to those compatible with the measured error syndrome. Our best strategy then is to classify the remaining, nondistinguishable errors according to their effect on the encoded logical information and undo the effects of the most probable equivalence class EE.

A hallmark of topological quantum error-correction codes [9-14] is the geometrical locality of these check operators: Physical qubits are placed on a lattice and check operators depend only on a few neighboring sites. The logical information, which is encoded globally in a subspace of all physical qubits, is preserved as long as we can successfully detect and correct local errors. If errors on the physical qubits occur with a probability p, the error threshold pc—a key figure of merit for any quantum code—is defined as the maximum error probability p, such that error classification is achievable. For error rates larger than pc, the error syndrome's complexity inhibits unambiguous error recovery. It is therefore of current interest to find codes where pc is large.

B. Error threshold as a phase transition

The process of error correction resembles a phase transition and, indeed, it is possible to connect error correction directly to an order-disorder phase transition in a suitable classical statistical-mechanical model [12,15,16]. One can derive a Hamiltonian HE of interacting Ising spins st, labeled by a Pauli error E that controls the sign of the couplings, such that the probability of each equivalence class EE is proportional to the partition function

p(EE) / ZE(p) := JV^^.

Equation (3) has to be interpreted as describing a random statistical model with quenched couplings and two parameters: the error probability p governing the fraction of negative interaction constants Ja 2 {±1}, and the inverse temperature / = 1/T. For low enough T and p the system orders into a ferromagnetic state (see Fig. 1). Along the Nishimori line [17] where Eq. (3) holds, the ordered (disordered) phase corresponds to the topological code being effective (ineffective). The intersection of the Nishimori line and the phase boundary identifies the error threshold pc.

The first topological codes studied were toric codes [9], still under intense investigation and scrutiny mainly due to their simplicity and elegance. To determine their error threshold, we show that toric codes under the depolarizing channel connect to the celebrated eight-vertex model (see Fig. 2) introduced by Sutherland [18], as well as Fan and Wu [19], and whose general solution by Baxter [20-22] stands as the culmination of a series of breakthroughs in the theory of phase transitions and critical phenomena.

The aforementioned mapping onto a statistical-mechanical model to compute the error tolerance of quantum codes was first applied to toric codes with bit-flip

Color codes (union jack)

0.00 0.05 0.10 0.15

Qubit error rate p

FIG. 1. Phase boundary estimated from Monte Carlo simulations for the estimation of the error threshold of the toric code, as well as two realizations of color codes (see text). The error threshold pc corresponds to the point at which the Nishimori lines intersects the phase boundary. Remarkably, the phase boundaries for all three codes agree within error bars. The stable ordered phase corresponds to the regime where quantum error correction is feasible.

FIG. 2. When computing the stability of the toric code to depolarization, the problem maps onto a classical statistical Ising model on two stacked square lattices. In addition to the standard two-body interactions for both top (a) and bottom (b) layers, the resulting Hamiltonian also includes four-body terms (c) that introduce correlations between the layers.

errors [15], connecting them to the random-bond Ising model. In general, for individual bit flips the error threshold is pc ~ 10.9%, and the same is true for phase flips alone. Therefore, under depolarizing noise and separately correcting bit flips and phase flips, the threshold is p'c = (3/2)pc ~ 16.4%. However, this result neglects correlations of bit flips and phase flips. We estimate the threshold under depolarizing noise for ideal error correction, such that, in particular, correlations are taken into account. We find pc = 18.9(3)%. Remarkably, the error threshold

increases significantly by taking correlation effects into account. They should thus not be neglected by recovery algorithms. A recent advance in this regard is the renor-malization approach of Duclos-Cianci and Poulin [23] where pc ' 16.4% was confirmed, still leaving room for further improvement [24]. Note also that pc is very close to the hashing bound p ' 0.1893 [26], which is also the case for uncorrelated bit-flip and phase-flip noise [15,27].

II. TOPOLOGICAL STABILIZER CODES A. Error correction in stabilizer codes

Both toric codes [9] and color codes [10] are topological stabilizer codes. A stabilizer code is described by a set of check operators S¿ in the Pauli group. That is, they are tensor products of Pauli operators and These

check operators S¿ are a commuting set of observables with eigenvalue ± 1 that generates an Abelian group S : = (S¿) that does not contain —1, called the stabilizer group. Encoded states | C) are those for which all check operators satisfy S¿|C) = +IC). If errors affect the state, typically they will change the value of the check operators, leaving a trace that can be used to recover the original state. Note that some errors are undetectable because they commute with all check operators and thus leave no trace.

We are interested in noisy channels of the form

P1 ! P2 = PoP + XP¿D¿PDy>

Po ! Pi = XP(E)EPoEt>

where the sum is over all Pauli group elements E, and p(E) denotes the probability for E to occur. Several different Pauli errors E have the same effect on the encoded state. Therefore, it is convenient to place them in equivalence classes E, such that E is equivalent to E' when EpEy = E'pE'y on an encoded state p or, equivalently, when EE' is proportional to a product of check operators. Therefore, the total probability for a given class of errors is given by

p(E) = £ p(SE).

One can choose a set of undetectable errors D,- and use them to label the error classes compatible with any given syndrome. Namely, if E is compatible with the syndrome, then the possible error classes are E itself and the classes DE.

The error-correction process starts with the measurement of the check operators S,-. Measuring each S,- yields an eigenvalue s,- = ±1. Only certain errors are compatible with these eigenvalues. In particular, E is compatible with the error syndrome if ES,- = SjSjE. Ideally, given a syndrome s = {s,-}, one can compute the relative probabilities P(E|s) of the different error classes E compatible with s. If Es is the class that maximizes this probability, the best guess is that this is the error that occurred and thus should be corrected. The net effect of such an ideal error correction is

where the success probability p0 and the probability for an effective error D,- are

:=XP(Es), Pi :=XP(DE). (7)

Note that in Eq. (7) the sum is over possible syndromes. Furthermore,

1 P(Es) , 1 ,

— — -:—r~ — 1, — — P0 — 1,

D P(s) D yo '

where D is the number of error classes per error syndrome. In practice, this ideal error correction might be too costly from a computational perspective. Therefore, approximations are needed.

B. Toric codes and color codes

Topological codes have two interesting features: First, they can be defined for different system sizes in such a way that check operators remain local—involving only a few neighboring qubits—and, at the same time, nontrivial un-detectable errors are global and thus involve a number of qubits that depend on the system size. Second, they exhibit an error threshold. For error rates below the threshold, the success probability [Eq. (7)] approaches 1 for increasing system size, whereas 1 — p0 decreases exponentially.

In toric codes [9], physical qubits are placed on the edges of a square lattice. Notice that, for each edge in the direct lattice, there is an edge in the dual lattice. Check operators Sf are attached to faces f, in either the direct or the dual lattices. Toric codes can thus be defined in two similar, but distinct ways: In the original definition by Kitaev [9], if f is a face in the direct (dual) lattice composed by the edges r, s, i, and w, then the corresponding check operator is Sf := ® ® ® (Sf : = ^r ® ® ^f ® The second definition, which is from Wen [28], does not distinguish between dual and direct faces. If f has a top edge r, a bottom edge s, and side edges t, w, then we take Sf := ^r ® ^f ® ® . Both definitions are equivalent up to a rotation of half of the qubits. However, for the depolarizing channel, Kitaev's definition is related to the alternating eight-vertex model and Wen's definition to the standard eight-vertex model.

In color codes [10], physical qubits are placed on the vertices of a trivalent lattice with three-colorable faces, such as, for example, the honeycomb lattice. There are two check operators Sf and Sf attached to each face f, taking the form Sf := <£)¿^f and Sf := <3>¿^f, respectively, with , running over the qubits on the vertices of f.

Because the computing capabilities of color codes depend on the underlying lattice where the qubits are placed, we study two different scenarios: the honeycomb lattice for its simplicity, and a lattice of octagons and

lattice are at the center of the plaquettes of the other. The Hamiltonian takes the explicit form

FIG. 3. For the hexagonal arrangement, there is a stabilizer operator Z6 for each of the hexagon plaquettes (a). In the mapping, these stabilizer operators translate to classical Ising spins, which are placed on the dual lattice [regular triangular lattice, (b)]. The square-octagonal setup (c) has wider computing capabilities because it allows for a larger class of quantum gates to be implemented. There are stabilizers Z4 (Z8) on the rectangles (octagons). The corresponding dual lattice in the mapping is the union-jack lattice (d).

squares that allows for the implementation of additional types of quantum gates. In the mapping onto a statistical-mechanical model to compute the error threshold, these two arrangements correspond to the triangular and union jack lattices, respectively (see Fig. 3).

III. RANDOM EIGHT-VERTEX ISING MODELS

To determine the error threshold, we show that topologi-cal codes under the depolarizing channel connect to certain random classical-spin models.

For the toric code, the error-correction process maps onto a statistical-mechanical interacting eight-vertex model [18-21]. Remarkably, this class of models exhibits critical exponents that depend on the coupling constants, Eq. (22) in Sec. IV B, thus challenging the very notion of universality. Eight-vertex models were originally formulated in the ''electric picture,'' where the degrees of freedom are electric dipoles placed at the bonds surrounding each vertex of a square lattice [22]; i.e., the number of independent dipole configurations per vertex is eight. In addition, a mapping to a ''magnetic picture'' was found by Wu [29], as well as Kadanoff and Wegner [30]: Consider two independent Ising systems, each on a square lattice, with classical spin variables st and taking on values ±1, and bonds Jj and J[', respectively. The lattices are stacked as shown in Fig. 2 such that the vertices (spin sites) of one

H = + JMSk4 +

This can be thought of as two interacting Ising models by means of a four-spin interaction (denoted by the symbol +) between original and dual lattices.

In fact, two types of eight-vertex models can be related to error correction in toric codes: the standard eight-vertex model, where Jj = J (Jj = J0) if a bond is a horizontal (vertical) link [see Figs. 2(a) and 2(b)], and the alternating eight-vertex model, where Jj = J (Jj = J0) if a bond belongs to the direct (dual) lattice. In both cases, we set the 4-spin interaction to J+ = J4, as depicted in Fig. 2(c). Thus, both types of eight-vertex models share the same lattice structure but differ in the pattern of coupling constants. This difference has fundamental consequences in the exact solvability of the model: While the standard eight-vertex model is exactly solved for arbitrary couplings, the alternating eight-vertex model is generally not exactly solvable. In fact, the latter corresponds to the Ashkin-Teller model [31]. Notice that, when J4 = 0, the eight-vertex model reduces to two decoupled Ising models, while for J4 ^ 0, the model has two critical temperatures.

The error threshold for correction in quantum codes corresponds to the critical line separating ordered from disordered phases. The ordered phase represents a situation in which quantum error correction can be performed with arbitrary precision. Determining the location of this critical line in eight-vertex models is facilitated by the existence of a self-duality symmetry in the partition function: a duality transformation relating a high-temperature eight-vertex model to a low-temperature one on the same lattice. Self-duality implies that the coupling constants for 2-spin interactions are isotropic, i.e., J = J0. Altogether, an iso-tropic, self-dual eight-vertex model has a critical line given by [22]

J0 = J;

-2pj4 = sinh(2pj),

with the restriction that J4 < J. The point in the plane (J4 = Jc, J = Jc ) at which the self-dual line ceases to be critical is given by

pJc = 4log(3) « 0.2746.

This is already a remarkable and encouraging result because it yields a critical point that is approximately 60% larger than in the standard square-lattice, two-dimensional Ising model. Note that the error threshold for bit-flip or phase-flip errors in the Kitaev model is computed via a mapping to the aforementioned two-dimensional Ising model. In that case, the critical point can be computed from the relationship sinh(2/Jc) = 1, i.e., /Jc = 0.4406. Recall that the critical exponents depend continuously on the value of J4.

FIG. 4. For topological color codes, qubits are arranged on trivalent lattices (hexagonal or square-octagonal). These codes are mapped to triangular lattices (triangular and union jack, respectively) with plaquette interactions (a),(b) on each layer, as well as six-body interactions correlating the two layers (c).

In this work we extend the standard eight-vertex model by adding quenched disorder to the couplings between the spins. Given that, for the eight-vertex model, /Jc = 0.2746 is smaller than for the square-lattice Ising model, we can expect to find a larger error threshold for the depolarizing channel than for bit-flip or phase-flip errors.

In addition to depolarizing errors in the toric code, where the problem maps onto Eq. (9), we also study color codes (see Fig. 3). In this case, the underlying statistical-mechanical model to study the error stability to depolarizing errors is defined on a triangular lattice. There are two Ising variables—si and si—per site. For convenience, we introduce an artificial third variable s0i0 = sis0i. The Hamiltonian is then given by

H2 = — X JkSiSjSk + Jksisjsk + J0jkSi0sj0sk0). (12)

hi j k>

Equation (12) is illustrated in Fig. 4 where the top (bottom) layer corresponds to the Sj (si) Ising variables, with the corresponding 3-spin interaction term as shown in Fig. 4(a) [4(b)]. The third term in the Hamiltonian with 6-spin interactions is represented by Fig. 4(c).

When Jjk = 0 in Eq. (12), we obtain two independent, triangular, three-body Ising models. Interestingly, this model can be mapped onto an eight-vertex model on a kagome lattice [32]. Therefore, the color-code Hamiltonian H2 in Eq. (12) can be thought of as an interacting eight-vertex model (or coupled eight-vertex models). In this work, we consider two different lattice geometries: triangular and union jack (see Fig. 3).

IV. MAPPING A. Spin models for depolarizing noise

The goal is to relate the stability of a topological stabilizer code to depolarizing noise to the ordered phase of a

suitably chosen classical spin model. However, here we consider the more general qubit channel shown in Eq. (1). This channel adds transparency to the mapping and reveals the differences between Kitaev's and Wen's versions of the toric code with respect to error correction. When Eq. (1) is applied to each qubit in a code, the net result is a channel of the form presented in Eq. (4). In particular, the probability for each Pauli error is

p(E) = (1 - p)n

w=x, y, Z

where n is the total number of qubits and Ew is the number of appearances of aw in the tensor product that forms E.

The classical spin Hamiltonian is constructed as follows:

(1) Attach a classical spin si to each check operator Si.

(2) Associate with each single-qubit Pauli operator a an interaction term Jasf sJT ■ ■■ sN such that the spins si correspond to the check operators Si affected by a, i.e., such that Sia = — aSi.

(3) Attach to each coupling a sign ra = ± 1 dictated by the Pauli error E labeling the Hamiltonian, through the conditions aE = rff£a.

The resulting Hamiltonian takes the general form

a 'as1s2

where the sum is over all Pauli operators a and there are

only three different couplings Ja since we set Jaw := J

with w = x, y, z, and k as the qubit label.

For the mapping, we require the interaction constants to be

Jw = - log

p^pypz

4P 6pw(1 - p)'

w = x, y, z.

This relates the error probability in Eq. (13) to the Boltzmann factor for the ordered state, {si = 1}, given the interactions generated by E:

p(E) ^ e-PH£({s; = 1}).

To recover Eq. (3), just notice that replacing the error E ! E0 = SjE is equivalent to considering a different spin configuration in the original Hamiltonian:

HsEfSi}) = He({(1 - 25;.■ )si}).

Finally, to complete the mapping, the label E in the Hamiltonian must describe quenched randomness. In particular, the coupling configuration dictated by E appears with probability p(E). Equivalently, this means that, for every qubit k, the probability p(r£, t^, t|) for each configuration of coupling signs is given by

p(+1, +1, +1) = 1 - p, p(+1,-1,-1) = px

p(-1, +1,-1) = py, p(-1,-1, +1) = pz.

In the case of the depolarizing channel,

= = = p/3 and Jx = Jy = Jz = J- (19)

The resulting model has two parameters, p and 0/, with 0/ = 1 /T. For low p and T, the model orders ferromag-netically and along the Nishimori line,

,-40/ =

p/3 (1 - p)'

this order is equivalent to the noise being correctable. Indeed, comparing error-class probabilities amounts to computing free-energy differences:

P(DE) _ ^(0) P(E) Ze(0)

= : e-0A;(0,E).

In topological codes, we expect the existence of an error threshold pc—or several for different error types, but we do not need such generality. When p < pc in the limit of large systems, the success probability is expected to approach unity, i.e., p0 ! 1. Thus, due to Eq. (8), along the Nishimori line the free-energy difference is asymptotically infinite, because for any real t, P(A"(0, E) > t)! 1. Similarly, when p > pc, the success probability is expected to become minimal (p0 ! 1 /D) and thus the free-energy difference converges in probability zero. Therefore, for any t > 0, we have P(|A"(0, E)| < t)! 1. This shows that the free-energy differences A" are order parameters, and pc is the critical value of p along the Nishimori line. In the models of interest here, these are domain-wall free energies.

B. Models for toric and color codes

Let us now study what the above mapping, Eq. (14), produces when applied to toric codes and to topological color codes.

For the toric code, the single-qubit operators aand az produce two-body interactions, because each bit flip (phase flip) affects the stabilizer operators on two neighboring dual [direct] faces. The ay operators, which combine correlated spin-flip and phase-flip errors, introduce four-body interactions (see Fig. 2). The result is an alternating eight-vertex model with coupling signs that are parametrized by a Pauli error E, namely,

He = -X(/*T+ SiSj + + JyT^SiSjSkS). (22)

The classical spin variables s (s[) live on the top (bottom) layer of two stacked, two-dimensional, Ising square lattices with interaction constants Jx (Jz) (see Fig. 2). The two layers are shifted by half a lattice spacing, and the third term in Eq. (22) describes the combined four-spin interaction at each of the crossings marked with the + symbol. Note that, in Eq. (22), there is one qubit per cross +. For Wen's toric code, one recovers the standard eight-vertex model. In either case, there is a global Z2 X Z2 symmetry because one can flip all spins in each lattice separately without affecting the total energy.

In the case of color codes, there is one spin per face. The a1 and az single-qubit operators produce three-body interactions in Eq. (14), whereas ay operators produce six-body interactions. The Hamiltonian is then given by Eq. (12) but with coupling signs parametrized by a Pauli error E, namely,

H = — / Tw „w„w„w

<i",j,k> w=x,y,z

with sfsysZ = 1. Therefore, we obtain two stacked, triangular lattices having three- and six-body interactions (see Fig. 4), with the six-body interactions introducing correlations between both layers. In this case, the global symmetry is Z2 X Z2 X Z2 X Z2. Indeed, the sites can be colored with three colors in such a way that each triangles has a site of each color. Thus one can flip all spins for two given colors in either of the two lattices separately without affecting the total energy.

For p = 0 in Eqs. (22) and (23), self-duality predicts a critical temperature of Tc = 1/0/c ' 3.641, a value that we confirm numerically in our Monte Carlo simulations.

V. MONTE CARLO SIMULATIONS

We investigate the classical statistical spin models acquired in the mapping, Eqs. (22) and (23), via large-scale classical Monte Carlo simulations using the parallel tempering Monte Carlo technique [33].

In the parallel tempering Monte Carlo method, several identical copies of the system at different temperatures are simulated. In addition to local simple Monte Carlo (Metropolis) spin updates [34], one performs global moves in which the temperatures of two neighboring copies are exchanged. It is important to select the position of the individual temperatures carefully such that the acceptance probabilities for the global moves are large enough [35] and each copy performs a random walk in temperature space. This choice, in turn, allows each copy to efficiently sample the rough energy landscape, therefore speeding up the simulation enormously.

Detecting the transition temperature Tc(p) for different fixed amounts of disorder allows us to pinpoint the phase boundary in the p-T phase diagram. The error threshold pc is then given by the intersection of the phase boundary with the Nishimori line.

A. Observables and simulation details

For the toric code, it is expedient to partition the lattice into two sublattices such that the only interconnection is given by the four-body interactions of the Hamiltonian in Eq. (22). The ground state of the pure system is realized when the spins of each sublattice are aligned (but the alignment may be different, as the sign would cancel out in both the two- and the four-spin terms). In this case, the sublattice magnetization is a good order parameter,

m =-V S;,

where P denotes one of the sublattices. Similarly, we note that each layer of the triangular lattice for color codes is tripartite, with spins aligned in each sublattice for all realizations of the pure system's ground state. Hence, we can define an order parameter analogous to Eq. (24) for a suitable sublattice P0. Note that particular caution is required when implementing the periodic boundary conditions to ensure that these distinct sublattices are well defined. We can now use the magnetization defined in Eq. (24) to construct the wave-vector-dependent magnetic susceptibility,

(k)-N; ((X )2 >t

where (• • -)T denotes a thermal average and R is the spatial location of the spin St. From Eq. (25), we construct the two-point finite-size correlation function [36],

2sin(kmin/2) V|> m (kmin )]

[*m (0)]a

where [• • -]av denotes an average over disorder, and kmin — (2^/L, 0; 0) is the smallest nonzero wave vector. Near the transition, is expected to scale as

&./L - 1[L^(T - Tc)],

TABLE I. Simulation parameters: L is the linear system size, Nsa is the number of disorder samples, teq — 2b is the number of equilibration sweeps (system size times number of single-spin Monte Carlo updates), Tmin (Tmax) is the lowest (highest) temperature, and NT is the number of temperatures used. For the toric code, we use L — {12, 16, 18, 24, 32, 36g, while, for color codes, we use L — {12, 15, 18, 24, 30, 36g, following the coloring constraints that the system size must be a multiple of 3.

p L Nsa b T ± min T ± max NT

0.00 12- -16 5000 18 3.500 4.000 42

0.00 18- 24 1000 19 3.500 4.000 42

0.00 30- 36 500 20 3.500 4.000 42

0.04- 0.05 12- 16 5000 20 3.200 3.800 42

0.04- 0.05 18- 24 1000 21 3.200 3.800 42

0.04- 0.05 30- 36 500 22 3.200 3.800 42

0.08- -0.12 12- 16 5000 20 2.700 3.500 42

0.08- -0.12 18- 24 1000 22 2.700 3.500 42

0.08- -0.12 30- 36 500 24 2.700 3.500 42

0.15 12- 16 5000 20 2.300 3.200 42

0.15 18- 24 1000 22 2.300 3.200 42

0.15 30- 36 500 24 2.300 3.200 42

0.17- -0.20 12- 16 5000 21 1.500 2.800 42

0.17- -0.20 18- 24 1000 23 1.500 2.800 42

0.17- 0.20 30- 36 500 25 1.500 2.800 42

where X is a dimensionless scaling function. Because at the transition temperature, T = Tc, the argument of Eq. (27) becomes zero (and hence independent of L), we expect lines of different system sizes to cross at this point. If, however, the lines do not meet, we know that no transition occurs in the studied temperature range.

In all simulations, equilibration is tested using a logarithmic binning of the data. This is done by verifying that all observables averaged over logarithmically increasing amounts of Monte Carlo time are, on average, time independent. Once the data for all observables agree for three logarithmic bins, we deem the Monte Carlo simulation for that system size to be in thermal equilibrium. The simulation parameters can be found in Table I.

B. Sample results

For the pure system (p = 0), there is a sharp transition visible directly in the sublattice magnetization. The transition temperature Tc pure ~ 3.64 coincides with the value found analytically. For larger amounts of disorder, a transition can still be located precisely by means of the crossings in the two-point finite-size correlation function [Eq. (26)] for different system sizes. Sample data for a disorder strength of p = 0.170 (i.e., meaning that on average 17% of the physical qubits are ''broken'') are shown in Fig. 5, indicating a transition temperature of Tc (p) = 2.14(2). The error bars are calculated using a bootstrap analysis of 500 samples. There are small finite-size effects which are addressed by analyzing the intersection T*(L, 2L) of pairs of system sizes. We estimate the limit value for L by means of a linear fit in a 1 /L

FIG. 5. Crossing of the correlation function for the toric code with a disorder rate of p = 0.170. The data exhibit a clear crossing at a transition temperature of Tc(p) ~ 2.14(2). Corrections to scaling are still minimal at this disorder rate, but they increase closer to the error threshold. Inset: Close-up of the area where the crossing occurs. The conservative estimate for the transition temperature is indicated by the vertical shaded area.

plot—This is our estimate for the best value in the physically relevant thermodynamic limit. For disorder rates approaching the error threshold, corrections to scaling increase and a careful finite-size scaling analysis has to be performed to determine Tc [37]. At p = 0.189, the lines touch only marginally such that both the scenario of a crossing as well as no transition are compatible within error bars. This gives rise to the large error bars in the phase diagram (Fig. 1). For error rates p > pc, the lines do not meet, indicating that there is no transition in the temperature range studied.

The crossing of the critical line Tc (p) with the Nishimori line [Eq. (20)] determines the error threshold to depolarization. Our (conservative) estimate is pc = 0.189(3). Our results are summarized in Fig. 1, which shows the estimated phase diagram.

VI. DUALITY METHOD

An alternative approach to estimating the critical value pc is to use the duality method [38], originally developed within the context of spin glasses.

The critical point of a statistical model expressed only by local interactions between degrees of freedom can be analyzed using the duality method under the assumption of a unique transition temperature. The partition function Z[A] is then given by the local Boltzmann factor A< = exp(0/cos^<), where < 2 {0, 1} is the difference between adjacent spins such as cos(^<) = ±1. We define the principal Boltzmann factor A0 as the case where all spins are parallel. The partition function has to be invariant under a Fourier transform, i.e., Z[A] = Z[A*], where A* is a dual principal Boltzmann factor (via a Fourier transformation). In that case, the critical point is determined via the equality A0 = A0. This implies that all the components of the local Boltzmann factors are equal for several self-dual models such as the standard Ising model. Although this self-duality does not work a priori for systems with quenched disorder in the general case, the method can be applied in a special subspace called the Nishimori line [38]. The results can be improved by considering extended local Boltzmann factors over a restricted

FIG. 6. Clusters used to estimate the error threshold for the depolarizing channel. The blue lines and triangles denote quenched random variables rZ, and the red lines and triangles correspond to rX. The central site is the spin variable summed over. The outer sites represent the spin variables fixed in the up direction.

range of interactions [38,39] (see Fig. 6, which illustrates the used clusters). Because the resulting statistical-mechanical Hamiltonians for both the toric code and the topological color codes are self-dual, we can apply this efficient technique to obtain estimates (up to systematic deviations that depend on the clusters used) of the error threshold.

A. Zeroth-order approximation

The effects of the depolarizing channel on topological codes can be expressed by a spin-glass model with the partition function [40]

I Vj'V. Vj

(r1, Tz) _

exp{0/rx + 0/rZ cos^<~

+ 0/rxrZ cos^0 cos^<~}.

Tj 2 {±1} and rZj 2 {±1} are quenched random variables chosen from the distribution

P(r*, rz.) / e^H+Tj+Tjj

This model has a gauge symmetry in the subspace J = Jp which corresponds to the Nishimori line.

To determine the multicritical point, we replicate the partition function to take into account the quenched ran-

domness of the variables rj and rZ,, i.e.,

ZJA] =

'(Y Y ПA(rjrj У"

LVi «A.- <j> 7 J

where the brackets denote a configurational average. The local Boltzmann factor is then given by

(rx.,rz.)

A Ij Ij

where k distinguishes the specific configuration (<, <~"). The n-binary Fourier transformation gives the dual Boltzmann factor AjJ k. It follows [38,39] that An,0 = AjJ 0 determines the critical point along the Nishimori line. Taking the leading term in n, we obtain the error threshold for the depolarizing channel of the toric code as pc = 0.189 ... under the conditions / = Jp and 3 exp(-40/) = p/(1 — p) for the Nishimori line. Because the local Boltzmann factors for the topological color codes on both the hexagonal and square-octagonal lattice are the same, we obtain the same estimate for the error threshold.

B. First-order approximation using finite clusters

To reduce systematic errors, we consider finite-size clusters with four bonds on each square lattice for the toric

code, six triangles taken from each triangular lattice for the color codes on the hexagonal lattice, and four triangles from each union-jack lattice for color codes on the square-octagonal lattice (see Fig. 6). We compute the principal Boltzmann factors on the clusters, i.e.,

A(1) = An,0

as well as its dual An(Q via an n-binary Fourier trans-

formation. As before, the critical point along the Nishimori line is determined via A^Q = An^. Taking the leading order in n, we obtain for the error thresholds

pc = 0.1888 ... (toric code),

pc = 0.1914... (color code, hexagonal),

ACKNOWLEDGMENTS

We would like to thank H. Nishimori and D. Poulin for useful discussions. M. A. M.-D. and H. B. thank the Spanish MICINN Grant No. FIS2009-10061, CAM research consortium QUITEMAD Grant No. S2009-ESP-1594, European Commission PICC: FP7 2007-2013, Grant No. 249958, UCM-BS Grant No. GICC-910758. Work at the Perimeter Institute is supported by Industry Canada and Ontario MRI. H. G. K. acknowledges support from the Swiss National Science Foundation (Grant No. PP002-114713) and the National Science Foundation (Grant No. DMR-1151387). M.O. acknowledges financial support from Grant-in-Aid for Young Scientists (B) No. 20740218 by MEXT and Kyoto University's GCOE Program Knowledge-Circulating Society. The authors acknowledge ETH Zurich for CPU time on the Brutus cluster and the Centro de Supercomputacion y Visualisation de Madrid (CeSViMa) for access to the Magerit-2 cluster.

pc = 0.1878 ... (color code, square octogonal). (36)

Although there are small variations in the estimates, the estimates are all of the order of approximately 19% and in agreement with our results from Monte Carlo simulations.

VII. RESULTS AND CONCLUSION

We have shown that the stability under depolarizing noise of toric codes can be related to the existence of a magnetic phase in a random eight-vertex model. Similarly, color codes turn out to be related to a class of "interacting" eight-vertex models. We analyze the models resulting from the mapping via both large-scale parallel tempering Monte Carlo simulations [16,37] and the duality method [38,39]. By determining Tc(p) for different error probabilities p, we are able to determine the phase boundary in the p-T plane (Fig. 1). Both approaches confirm the existence of a stable ordered phase and by locating the intersection of the phase boundary with the Nishimori line, we compute, in a nonperturbative way, the disturbing effects of noise on topological codes. The external noise considered in this work is more realistic than in previous studies because it applies to bit-flip errors, phase-flip errors, and, more importantly, a nontrivial combination thereof.

The error threshold to depolarization errors for different classes of topological codes studied is approximately 19%, which is larger than the threshold for noncorrelated errors. This is very encouraging and shows that topological codes are more resilient to depolarization effects than previously thought. The profound relationship between complex statistical-mechanical models, such as the eight-vertex model, and topological quantum error correction promises to deliver a plethora of new paradigms to be studied in both fields in coming years.

[1] P. W. Shor and J. Preskill, Simple Proof of Security of the BB84 Quantum Key Distribution Protocol, Phys. Rev. Lett. 85, 441 (2000).

[2] B. Kraus, N. Gisin, and R. Renner, Lower and Upper Bounds on the Secret-Key Rate for Quantum Key Distribution Protocols Using One-Way Classical Communication, Phys. Rev. Lett. 95, 080501 (2005).

[3] C. H. Bennett, G. Brassard, S. Popescu, B. Schumacher, J. A. Smolin, and W. K. Wootters, Purification of Noisy Entanglement and Faithful Teleportation via Noisy Channels, Phys. Rev. Lett. 76, 722 (1996).

[4] G. Bowen and S. Bose, Teleportation as a Depolarizing Quantum Channel, Relative Entropy, and Classical Capacity, Phys. Rev. Lett. 87, 267901 (2001).

[5] A. Shaham and H. S. Eisenberg, Realizing Ccontrollable Depolarization in Photonic Quantum-Information Channels, Phys. Rev. A 83, 022303 (2011).

[6] A. Chiuri, V. Rosati, G. Vallone, S. Padua, H. Imai, S. Giacomini, C. Macchiavello, and P. Mataloni, Experimental Realization of Optimal Noise Estimation for a General Pauli Channel, Phys. Rev. Lett. 107, 253602 (2011).

[7] P. W. Shor, Scheme for Reducing Decoherence in Quantum Computer Memory, Phys. Rev. A 52, R2493 (1995).

[8] A. M. Steane, Error Correcting Codes in Quantum Theory, Phys. Rev. Lett. 77, 793 (1996).

[9] A. Y. Kitaev, Fault-Tolerant Quantum Computation by Anyons, Ann. Phys. (N.Y.) 303, 2 (2003).

[10] H. Bombin and M. A. Martin-Delgado, Topological Quantum Distillation, Phys. Rev. Lett. 97, 180501 (2006).

[11] H. Bombin and M. A. Martin-Delgado, Exact Topological Quantum Order in D = 3 and Beyond: Branyons and Brane-Net Condensates, Phys. Rev. B 75, 075103 (2007).

[12] H. Bombin, Topological Subsystem Codes, Phys. Rev. A 81, 032301 (2010).

[13] S. Bravyi, B.M. Terhal, and B. Leemhuis, Majorana Fermion Codes, New J. Phys. 12, 083039 (2010).

[14] J. Haah, Local Stabilizer Codes in Three Dimensions without String Logical Operators, Phys. Rev. A 83, 042330 (2011).

[15] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, Topological Quantum Memory, J. Math. Phys. (N.Y.) 43, 4452 (2002).

[16] H. G. Katzgraber, H. Bombin, and M. A. Martin-Delgado, Error Threshold for Color Codes and Random 3-Body Ising Models, Phys. Rev. Lett. 103, 090501 (2009).

[17] H. Nishimori, Internal Energy, Specific Heat and Correlation Function of the Bond-Random Ising Model, Prog. Theor. Phys. 66, 1169 (1981).

[18] B. Sutherland, Two-Dimensional Hydrogen Bonded Crystals without the Ice Rule, J. Math. Phys. (N.Y.) 11, 3183 (1970).

[19] C. Fan and F. Y. Wu, General Lattice Model of Phase Transitions, Phys. Rev. B 2, 723 (1970).

[20] R. J. Baxter, Eight-Vertex Model in Lattice Statistics, Phys. Rev. Lett. 26, 832 (1971).

[21] R. J. Baxter, Partition Function of the Eight-Vertex Lattice Model, Ann. Phys. (N.Y.) 70, 193 (1972).

[22] R. Baxter, Exactly Solved Models in Statistical Mechanics (Academic Press, London, 1982).

[23] G. Duclos-Cianci and D. Poulin, Fast Decoders for Topological Quantum Codes, Phys. Rev. Lett. 104, 050504 (2010).

[24] After completion of this work, we learned of a new algorithm [25] based on parallel tempering Monte Carlo simulation that attains an error threshold of pc ' 18.5%.

[25] J. R. Wootton and D. Loss, High Threshold Error Correction for the Surface Code, arXiv:quant-phys/ 1202.4316.

[26] C. H. Bennett, D. P. Divincenzo, J. A. Smolin, and W. K. Wootters, Mixed-State Entanglement and Quantum Error Correction, Phys. Rev. A 54, 3824 (1996).

[27] B. Rothlisberger, J. R. Wootton, R. M. Heath, J. K. Pachos, and D. Loss, Incoherent Dynamics in the Toric Code Subject to Disorder, Phys. Rev. A 85, 022313 (2012).

[28] X.-G. Wen, Quantum Orders in an Exact Soluble Model, Phys. Rev. Lett. 90, 016803 (2003).

[29] F. W. Wu, Ising Model with Four-Spin Interactions, Phys. Rev. B 4, 2312 (1971).

[30] L. P. Kadanoff and F. J. Wegner, Some Critical Properties of the Eight-Vertex Model, Phys. Rev. B 4, 3989 (1971).

[31] J. Ashkin and E. Teller, Statistics of Two-Dimensional Lattices with Four Components, Phys. Rev. 64, 178 (1943).

[32] R. J. Baxter, Solvable Eight-Vertex Model on an Arbitrary Planar Lattice, Phil. Trans. R. Soc. A 289, 315 (1978).

[33] K. Hukushima and K. Nemoto, Exchange Monte Carlo Method and Application to Spin Glass Simulations, J. Phys. Soc. Jpn. 65, 1604 (1996).

[34] M.E.J. Newman and G.T. Barkema, Monte Carlo Methods in Statistical Physics (Oxford University Press, New York, 1999).

[35] H. G. Katzgraber, M. Körner, and A. P. Young, Universality in Three-Dimensional Ising Spin Glasses: A Monte Carlo Study, Phys. Rev. B 73, 224432 (2006).

[36] M. Palassini and S. Caracciolo, Universal Finite-Size Scaling Functions in the 3D Ising Spin Glass, Phys. Rev. Lett. 82, 5128 (1999).

[37] R. S. Andrist, H. G. Katzgraber, H. Bombin, and M. A. Martin-Delgado, Tricolored Lattice Gauge Theory with Randomness: Fault-Tolerance in Topological Color Codes, New J. Phys. 13, 083006 (2011).

[38] M. Ohzeki, Precise Locations of Multicritical Points for Spin Glasses on Regular Lattices, Phys. Rev. E 79, 021129 (2009).

[39] M. Ohzeki, Accuracy Thresholds of Topological Color Codes on the Hexagonal and Square-Octagonal Lattices, Phys. Rev. E 80, 011141 (2009).

[40] H. Nishimori and K. Nemoto, Duality and Multicritical Point of Two-Dimensional Spin Glasses, J. Phys. Soc. Jpn. 71, 1198 (2002).