Tricolored Lattice Gauge Theory with Randomness: Fault-Tolerance in Topological Color Codes

We compute the error threshold of color codes, a class of topological quantum codes that allow a direct implementation of quantum Clifford gates, when both qubit and measurement errors are present. By mapping the problem onto a statistical-mechanical three-dimensional disordered Ising lattice gauge theory, we estimate via large-scale Monte Carlo simulations that color codes are stable against 4.5(2)% errors. Furthermore, by evaluating the skewness of the Wilson loop distributions, we introduce a very sensitive probe to locate first-order phase transitions in lattice gauge theories.


Introduction
The study of quantum error correction in topological stabilizer codes [1] burgeoned a magnificent synergy between quantum computation and statistical mechanical systems with disorder [2]. Quantum error correction [3,4] is a method that preserves quantum information from decoherence by actively detecting and counteracting errors using redundancy. In topological codes, these error-detecting operations are local and one can relate the stability of a topological quantum memory to an ordered phase in a classical statistical model [2]. The error threshold-a figure of merit that describes the stability of a system against local errors and up to which error correction is feasible-can be identified with the critical point at which the ordered phase of the underlying statistical mechanical system is lost due to the influence of quenched disorder.
In most studies, it is assumed that the error-correction procedure is free of errors. However, error correction is achieved by means of applying a set of quantum gates and this procedure can be flawed as well, leading to the concept of fault tolerance [5][6][7]. In a topological quantum memory, the information is encoded in global degrees of freedom [2] and preserved by repeatedly performing local measurements to keep track of and correct for errors. To understand the error resilience, it is convenient to adopt a phenomenological description of errors [2]. Errors are assumed to be nonsystematic and uncorrelated both in space and in time. Therefore, the error modeling process is parameterized in terms of two error rates: the qubit error rate p and the measurement error rate q. At each complete step in the syndrome measurement process, each physical qubit in the memory can suffer an error with an independent probability p and each measurement outcome can be incorrect with probability q. For the case of toric codes [1], which were the first example of topological codes, error correction (q = 0) can be studied via the twodimensional (2D) random-bond Ising model [8,9]. Including faulty measurements (q > 0), the error process is mapped onto a 3D random-plaquette gauge model [2,10].
More recently [11,12], other approaches to topological quantum error correction have been introduced. For example, topological color codes (TCCs) are a class of topological codes that allow the transversal implementation of the Clifford group of quantum gates [11,12]. Remarkably, these enhanced computational capabilities for quantum distillation, teleportation, dense coding, etc are possible while preserving a high error threshold [13] (in comparison to the Kitaev model [1]). This result has been obtained numerically after mapping the error correction (q = 0) for TCCs onto a random three-body Ising model on a triangular lattice [13], later confirmed by other numerical methods [14][15][16]. It is, however, unclear whether quantum computations on TCCs can be performed reliably in the presence of faulty measurements (q = 0). Here we address this problem by simulating a complex disordered lattice gauge theory (LGT [17])-a task that requires the introduction of novel tools to probe for an ordered phase. Our main results are: (1) A complete study of error correction in TCCs. We find that, for equal error qubit flip and measurement error rates p = q, the error threshold is p c = 4.8(2)%. (2) A novel Abelian LGT with gauge group Z 2 × Z 2 and a peculiar lattice and gauge structure that departs from the standard formulations of Wegner [18] and Wilson [19]; see figure 1. We refer to it as a tricolored LGT. Its structure reflects the error history in color codes rather than the discretization of a continuous gauge theory. (3) A novel approach to pinpoint first-order phase transitions in LGTs with disorder using the skewness of the average over Wilson loop operators.

Error correction with faulty measurements
To construct a TCC, we consider a hexagonal arrangement of the physical qubits on a 2D surface, such that the dual lattice has a regular triangular form with three-colorable vertices (i.e. the vertices can be labeled with three colors such that adjacent vertices have different colors). Note that each triangle of the dual lattice T (see figure 1) corresponds to a physical qubit in the initial arrangement. We now attach to each vertex v two Pauli operators, which act on the physical qubits associated with the six triangles adjacent to v. These are called check operators because encoded states |ψ satisfy X v |ψ = Z v |ψ = |ψ . Such states exist because the group generated by check operators, called a stabilizer group, does not contain −1, so that, in particular, check operators commute with each other. The dimension of the encoded subspace depends only on the topology of the surface where the code lives, something that is true in part due to the colorability properties of the lattice. For example, a regular lattice with periodic boundary conditions has the topology of a torus and encodes two logical qubits [11]. As in any topological code, the main feature is that local operators cannot distinguish different encoded states.
To keep track of errors, check operators are measured repeatedly over time, allowing for the detection of local inconsistencies with the code. Error correction amounts to guessing the correct error history E among those that are compatible with the measurement outcomes. Indeed, many error histories have an equivalent effect, and thus the ideal strategy is to compute which equivalence classĒ happened with the highest probability P(Ē). Therefore, error correction 4 is highly successful when for typical errors there is a class that dominates over the others. Color codes, being topological, have an error threshold below which the success probability approaches unity in the limit of large codes.
Because color codes are Calderbank-Shor-Steane (CSS) codes [20,21], X (bit-flip) and Z (phase-flip) errors can be corrected separately. In doing so one ignores any correlations between bit-flip and phase-flip errors, something that can only make error correction worse. Then, for simplicity, we assume an error model where bit-flip and phase-flip errors are uncorrelated. Before each round of measurements, each qubit is subject to the channel where ρ represents the state of a single qubit and p is both the probability for a bit-flip error to occur and the probability for a phase-flip error to occur. As for the measurement outcomes of check operators, each of them is wrong with independent probability q.
We focus on bit-flip errors, which change the eigenvalue of all adjacent Z v check operators; Z errors are analogous. Identifying time with the vertical direction, we form a lattice of stacked triangle layers, one per measurement round, with vertical connections (see figure 1, left). Error histories are represented as collections of triangles representing bit-flip errors and vertical links that represent measurement errors. If an error history E is composed of a total of b bit-flip errors and m wrong measurement outcomes, its probability is where N is the total number of qubits and M the total number of check operators. Following [2], we classify error histories E in equivalence classesĒ. Whenever two error histories E and E belong to the same class and an error E has happened, the error correction is still successful even if we wrongly correct according to E . The equivalence relation can be constructed from two different kinds of elementary equivalences. Firstly, if S is the set of triangles meeting at a given vertex v ( figure 1 (a)) and the histories E and E differ exactly on the elements of S, then they are equivalent because the local operator X v has no effect on encoded states. In particular, if P is a Pauli operator representing an error history-including X and Z errors-and ρ an encoded state, then We say that these elementary equivalences are of type I. Secondly, suppose that S is now the set including two triangles, one on top of the other, and the vertical links connecting them ( figure 1(b)). Then if the histories E and E differ exactly on the elements of S, they are equivalent because two unnoticed subsequent bit-flips on the same qubit are irrelevant. We say that these elementary equivalences are of type II.

Tricolored lattice gauge theory with randomness
As in [2], the first goal is to express the probabilities of error classes in terms of the partition function of a suitable statistical model. Our goal is to obtain 5 for a suitably parameterized Hamiltonian H E . We start by attaching classical spins σ = ±1 to the elementary equivalences between error histories discussed above, both of types I and II. This produces a new lattice (see figure 1 (right)) with honeycomb layers H sandwiched between triangular layers T . Spins on the T -layers (H-layers) correspond to equivalences of type I (type II) while negative interaction constants correspond to errors. There are six-body interactions that stem from measurement errors (vertical links in the previous lattice) that involve the vertices of a given hexagon at an H-layer; see figure 1(d). There are also five-body interactions due to bit-flip errors and thus related to triangles in T -layers: they involve the three vertices of the triangle plus the two H-vertices directly above and below the triangle; see figure 1(c). We use the symbol · to denote such sets of five vertices. The Hamiltonian of interest is where J > 0 and K > 0, γ · , γ = ±1 and j [k] runs over the vertices of each · [ ]. Given an error history E, let γ E be such that γ E · = −1 when the triangle in · belongs to E, and similarly for hexagons and their dual vertical links. Also, let σ 0 represent the state with all spins in the +1 state. Then if the inverse temperature β = 1/T as well as J and K satisfy the conditions Moreover, if two spin configurations σ 1 and σ 2 only differ by the sign of a single spin, and two error histories E 1 and E 2 are equivalent up to the corresponding elementary equivalence, then Putting together equations (7) and (9), we obtain the desired result Following [2], the next step is to consider a model with a Hamiltonian (equation (6)) where the parameters γ are quenched variables, with p and q dictating the probability distribution of the γ 's such that P(γ E ) := P(E). That is, each · [ ] has negative sign with probability p [q]. The resulting system has four parameters, β J , β K , p and q. The connection between error correction in color codes and this statistical model happens along the sheet described by the Nishimori conditions (equations 8). In that sheet, order in the statistical model corresponds to errors being correctable [2]. For the sake of simplicity, we set p = q, assuming the same fault rate for qubit and measurement errors. In that case it is enough to consider a statistical model with K = J , so that it has only two parameters, namely β J and p. The critical p along the Nishimori line then gives the fault-tolerance threshold p c for color codes.

Gauge symmetry
The Hamiltonian in equation (6) has a local symmetry per hexagon : flipping its spins and those above and below its center (figure 1(e)) leaves the energy invariant. This is a Z 2 × Z 2 gauge 6 symmetry, as the construction of suitable Wilson loops demonstrates. Let us label H-hexagons according to the coloring of the T -vertices. Wilson loops also come in three colors. For example, a blue horizontal loop is the product of green and blue hexagonal terms ( figure 1(f)). Since there are two independent Wilson loops for a given surface with a single boundary, there are four possible flux values going through it. Moreover, the total flux through two regions is trivial if both have the same flux, showing that the gauge is Z 2 × Z 2 . Indeed, excitations can be regarded as colored fluxes. Firstly, depict an excited hexagon as a vertical flux of the corresponding color. Secondly, depict an excited triangle as the three merging fluxes (one of each color). Using this convention, excitations take the form of closed colored strings with branching points where three different colors meet ( figure 1(g)).

Simulation details
The study of the partition function (equation (6)) constitutes a considerable numerical challenge.
To compute the error threshold including bit-flip and measurement errors, we need to compute the p -q -T c phase diagram of the model. The special case of equal error rates ( p = q) studied here provides useful guidance; however, 23.4 CPU years (6.4 × 10 15 operations on Brutus) were needed to obtain the results in figure 5. Because the model is an LGT, the natural (gauge invariant) order parameter is the Wilson loop. Owing to the presence of disorder and the complexity of the lattice, a clean scaling analysis using the area/perimeter laws [10] is imprecise because large Wilson loops show strong corrections to scaling. Instead, we investigate in detail the smallest horizontal loops in the system, which correspond to single hexagon plaquettes , and record their average value over all loops in the system, i.e.
suitably averaged over Monte Carlo time and disorder. Without disorder ( p = 0) this order parameter shows a clear signal of a transition. However, detection of the transition temperature T c becomes increasingly difficult when a finite error probability ( p > 0) is introduced. Therefore, we also study the whole distribution f (w). For a first-order phase transition, we expect a characteristic double-peak structure to emerge near the transition. To better pinpoint the transition, we introduce a derived order parameter ζ related to the skewness of the distribution wherew = w − [ w T ] av , · · · T denotes a thermal average and [· · ·] av represents an average over the disorder. The skewness of a symmetric function is zero and becomes positive (negative) when the function is slanted to the left (right). Therefore, we can identify the point where ζ crosses the horizontal axis with an equally weighted double-peak structure in f (w), i.e. ζ (T = T c ) = 0. We have compared these results to the peaks in the specific heat as well as a Maxwell construction and find perfect agreement. To compute the p -T c phase diagram, we fix the value of the error rate p and vary the temperature T until we detect a Higgs-to-confining Table 1. Simulation parameters: L is the layer size, M is the number of layers, N sa is the number of disorder samples, t eq = 2 b is the number of equilibration sweeps, T min (T max ) is the lowest (highest) temperature and N T is the number of temperatures used. transition for T = T c . T c ( p) is then the critical line separating both phases. The error threshold p c is given by the crossing point between T c ( p) and the Nishimori line [22]. The term 'Higgs-to-confining' refers to the transition from one gapped phase to another in an LGT. These gapped phases are characterized by non-vanishing closed string correlators, the Wilson loops, and can be distinguished by investigating the perimeter and area laws, respectively. The Higgs phase corresponds to a perimeter law for the Wilson loop correlator, regardless of whether an explicit Higgs field is present in the theory. The rationale is that static matter sources, when introduced in the Higgs phase, will become unbounded, i.e. deconfined, as opposed to the confinement that occurs in the phase with an area law. Therefore, here the Higgs phase refers to the (magnetically) ordered phase, whereas the confined phase corresponds to the disordered phase found for larger values of p and T .
It has been previously shown [2] that the error threshold is given by the crossing point of the Nishimori line with this phase boundary. We recall that the Nishimori line T = T ( p) is the locus describing a quantum computer in the presence of external noise, while the rest of the phase diagram is merely introduced as an auxiliary tool to locate the multicritical point. Note that on the Nishimori line, the effects of thermal fluctuations and quenched randomness are in balance. For weak disorder p and low temperatures T , the system is in a magnetically ordered Higgs phase. In terms of the color code, this indicates that all likely error histories for a given error syndrome are topologically equivalent and error recovery is achievable. However, at a critical disorder value p c (and a corresponding temperature determined by Nishimori's relation), magnetic flux tubes condense and the system enters the magnetically disordered confinement phase. In this case, magnetic disorder means that the error syndrome cannot point to likely error patterns belonging to a unique topological class; therefore, the topologically encoded information is vulnerable to damage.
For the simulations, we study systems of size L 2 × M with M along the vertical direction. Periodic boundary conditions are applied to reduce finite-size effects. The colorability and stacking of the layers require that L [M] is a multiple of 3 [2]. For the system to be as isotropic as possible, we choose parameter pairs such that M ∈ {L , L − 1}. Because the numerical complexity of the system increases considerably when p > 0, we use the parallel tempering Monte Carlo technique [23]. Simulation parameters are listed in table 1. Equilibration is tested by a logarithmic binning of the data. Once the last three bins agree within statistical error bars, the system is deemed to be in thermal equilibrium.   Figure 2 shows the average Wilson loop value as a function of temperature for the pure system ( p = 0) and for an error rate of p = 0.03 (inset). For the system without disorder, the transition is clearly visible (note that extrapolating to the thermodynamic limit shows that for p = 0 the transition seems to also be first order; see also figure 4). However, when p > 0 it is difficult to determine the location of the transition. An alternative view is provided by the histogram of Wilson loop values f (w) for different temperatures ( figure 3). Below T c (left panel), one can observe the development of a shoulder that eventually becomes a second peak. The two peaks have a symmetric weight distribution at the transition temperature (center panel). Above T c the initial peak starts to disappear (right panel) and the distribution slants to the right. This property is mirrored by the skewness of the distribution (figure 4). ζ (w) has a positive (negative) peak where the second (first) peak in f (w) develops (disappears), with a zero where the weight distribution is symmetric, i.e. at T = T c . We have compared our results using the skewness to the peak position of the specific heat, as well as a Maxwell construction. Not only does the skewness deliver the most precise results, for large disorder ( p 0.03) it is the only method that reliably shows a sign of a transition, if present. We estimate T c in the thermodynamic limit by plotting the size-dependent crossing point T * c (N ) against 1/N and applying a linear fit (inset). The full phase diagram is shown in figure 5; the crossing between T c ( p) and the Nishimori line (thin blue line) determines the error threshold p c = 0.048 (2). Finally, note that both vertical and horizontal Wilson loops give T c values that agree within error bars.

Conclusions
By mapping TCCs with both qubit and measurement errors onto a classical statistical mechanical model with disorder, we have obtained an Abelian LGT with a peculiar lattice and gauge structure. To date, toric codes were the only codes whose complete fault-tolerance properties had been studied. Using Monte Carlo simulations, we estimate the error threshold for TCCs with both bit-flip and measurement errors to be p c = 4.8(2)% (to be compared with p c ≈ 3% for the toric code [10]). To obtain this result we introduce a new approach that uses the skewness of the Wilson loop distribution to pinpoint first-order phase transitions in lattice gauge theories with disorder.