Topological quantum error correction in the Kitaev honeycomb model

TL;DR

This paper analyzes the topological quantum error correction capabilities of the Kitaev honeycomb model, providing rigorous bounds and numerical insights, and discovering a potential lifetime scaling advantage over the toric code in certain regimes.

Contribution

It offers the first explicit exponential bounds on degeneracy, indistinguishability, and correctability for the honeycomb model, and explores its error correction performance near the toric code fixed point.

Findings

Exponential bounds on code properties are tighter than general topological phase bounds.

Thermal noise has minimal impact on error correction in studied regimes.

Low-temperature, finite-size effects may improve lifetime scaling compared to the toric code.

Abstract

The Kitaev honeycomb model is an approximate topological quantum error correcting code in the same phase as the toric code, but requiring only a 2-body Hamiltonian. As a frustrated spin model, it is well outside the commuting models of topological quantum codes that are typically studied, but its exact solubility makes it more amenable to analysis of effects arising in this noncommutative setting than a generic topologically ordered Hamiltonian. Here we study quantum error correction in the honeycomb model using both analytic and numerical techniques. We first prove explicit exponential bounds on the approximate degeneracy, local indistinguishability, and correctability of the code space. These bounds are tighter than can be achieved using known general properties of topological phases. Our proofs are specialized to the honeycomb model, but some of the methods may nonetheless be of broader interest. Following this, we numerically study noise caused by thermalization processes in the perturbative regime close to the toric code renormalization group fixed point. The appearance of non-topological excitations in this setting has no significant effect on the error correction properties of the honeycomb model in the regimes we study. Although the behavior of this model is found to be qualitatively similar to that of the standard toric code in most regimes, we find numerical evidence of an interesting effect in the low-temperature, finite-size regime where a preferred lattice direction emerges and anyon diffusion is geometrically constrained. We expect this effect to yield an improvement in the scaling of the lifetime with system size as compared to the standard toric code.

Equations

$$\Delta E\le c_1{N}^2{\mathrm{e}}^{-aN},$$

$$\left|{⟨O_A⟩}_{{\psi }_1}-{⟨O_A⟩}_{{\psi }_2}\right|\le c_3{\mathrm{e}}^{-aN/4}\Vert O_A\Vert ,$$

$$ϵ(N)\le c_5{\mathrm{e}}^{-aN/16},$$

$$H=-\sum _{\mathbf{q}}(J_x{K}_{\mathbf{q}}^x+J_y{K}_{\mathbf{q}}^y+J_z{K}_{\mathbf{q}}^z)$$

$$K_{\mathbf{q}}^x$$

$$K_{\mathbf{q}}^y$$

$$K_{\mathbf{q}}^z$$

$$W_{\mathbf{q}}$$

$$={\sigma }_{\mathbf{q},\circ }^z{\sigma }_{\mathbf{q}+{\mathbf{n}}_y,\bullet }^x{\sigma }_{\mathbf{q}+{\mathbf{n}}_y,\circ }^y{\sigma }_{\mathbf{q}+{\mathbf{n}}_x+{\mathbf{n}}_y,\bullet }^z{\sigma }_{\mathbf{q}+{\mathbf{n}}_x,\circ }^x{\sigma }_{\mathbf{q}+{\mathbf{n}}_x,\bullet }^y$$

$$L_x$$

$$L_y$$

$$J_z>J_x\ge J_y>0\text{and}{J}_z>J_x+J_y,$$

$$K_{\mathbf{q}}^x$$

$$K_{\mathbf{q}}^y$$

$$K_{\mathbf{q}}^z$$

$$H_{\mathrm{vf}}=\sum _{\mathbf{k}}\left[{\xi }_{\mathbf{k}}{c}_{\mathbf{k}}^{†}{c}_{\mathbf{k}}+\frac{1}{2}({\Delta }_{\mathbf{k}}{c}_{\mathbf{k}}^{†}{c}_{-\mathbf{k}}^{†}+{\Delta }_{\mathbf{k}}^{\ast }{c}_{-\mathbf{k}}{c}_{\mathbf{k}})\right]-N^2{J}_z,$$

$$c_{\mathbf{k}}$$

$${\xi }_{\mathbf{k}}$$

$${\Delta }_{\mathbf{k}}$$

$$H_{\mathrm{vf}}=\sum _{\mathbf{k}}{E}_{\mathbf{k}}\left({\gamma }_{\mathbf{k}}^{†}{\gamma }_{\mathbf{k}}-\frac{1}{2}\right)$$

$${\gamma }_{\mathbf{k}}=u_{\mathbf{k}}{c}_{\mathbf{k}}-v_{\mathbf{k}}{c}_{-\mathbf{k}}^{†},$$

$$u_{\mathbf{k}}$$

$$v_{\mathbf{k}}$$

$$E_{\ell }=\sum _{\mathbf{k}}\left(-\frac{1}{2}{E}_{\mathbf{k}}\right),$$

$$E_{\mathbf{k}}$$

$$k_{\alpha }$$

$$k_x\stackrel{N\to \infty }{\to }x\in [0,2\pi ),$$

$$\frac{E_{\ell }}{N^2}\stackrel{N\to \infty }{\to }\frac{1}{4{\pi }^2}{\int }_0^{2\pi }{\int }_0^{2\pi }-\frac{1}{2}E(x,y)\mathrm{d}x\mathrm{d}y.$$

$$V=-\sum _{\mathbf{q}}(J_x{K}_{\mathbf{q}}^x+J_y{K}_{\mathbf{q}}^y)$$

$$H_{\mathrm{eff}}=U^{†}HU$$

$$H_{\mathrm{eff}}=E_0+\mu H_0-\sum _{\{{\mathbf{q}}_1,\dots ,{\mathbf{q}}_n\}}{C}_{\{{\mathbf{q}}_1,\dots ,{\mathbf{q}}_n\}}{W}_{{\mathbf{q}}_1}\cdots W_{{\mathbf{q}}_n}-\sum _{\{{\mathbf{q}}_1^{\prime },\dots ,{\mathbf{q}}_n^{\prime }\}}{D}_{\{{\mathbf{q}}_1^{\prime },\dots ,{\mathbf{q}}_n^{\prime }\}}{R}_{\{{\mathbf{q}}_1^{\prime },\dots ,{\mathbf{q}}_n^{\prime }\}}$$

$${\Pi }_0{H}_{\mathrm{eff}}{\Pi }_0=E_0^{\prime }-\sum _{\{{\mathbf{q}}_1,\dots ,{\mathbf{q}}_n\}}{C}_{\{{\mathbf{q}}_1,\dots ,{\mathbf{q}}_n\}}{\Pi }_0{W}_{{\mathbf{q}}_1}\cdots W_{{\mathbf{q}}_n}{\Pi }_0.$$

$${\Pi }_0{H}_{\mathrm{eff}}^{(4)}{\Pi }_0=E_0^{\prime }-\sum _{\mathbf{q}}{C}_{\mathbf{q}}{\Pi }_0{W}_{\mathbf{q}}{\Pi }_0.$$

$$H(s)$$

$$\Vert [{\mathrm{e}}^{iH(s)t}{O}_A{\mathrm{e}}^{-iH(s)t},O_B]\Vert \le C\Vert O_A\Vert \Vert O_B\Vert \min (|A|,|B|){\mathrm{e}}^{-\alpha \left(d(A,B)-v|t|\right)}$$

Full text

I Introduction

Topologically ordered systems are a promising avenue to implement quantum information storage in a physical system Kitaev (2003); Dennis et al. (2002). These systems have degenerate ground spaces in which information can be encoded such that it is robust to local perturbations Bravyi et al. (2010); Bravyi and Hastings (2011); Michalakis and Zwolak (2013). These kinds of systems can also be used to perform topologically protected quantum computation Nayak et al. (2008). The prototypical example of a topological spin model is the toric code Kitaev (2003). This model involves local interactions in 2D, and gives rise to abelian anyonic quasiparticles. The toric code has many desirable features, such as exact degeneracy of the Hamitonian eigenspaces for all system sizes, harmonic spectrum, flat dispersion relation, and ground states that are renormalization group (RG) fixed points Levin and Wen (2005); Aguado and Vidal (2008). This wealth of special structure makes the analysis Dennis et al. (2002); Alicki et al. (2009); Viyuela et al. (2012); Freeman et al. (2014); Jouzdani et al. (2014) and simulation Brown et al. (2016); Darmawan and Poulin (2017) of noise and error correction protocols tractable, at least for certain types of noise. However, the toric code is experimentally challenging to implement due to the fact that it involves 4-body interactions. Additionally, the special structure that enables such detailed analysis might obscure some features of error correction and noise that may be present in more general models.

The Kitaev honeycomb model Kitaev (2006) is a topologically ordered system involving only 2-body interactions, and is considered an experimentally plausible model to implement in many physical systems Duan et al. (2003); Micheli et al. (2006); Zhang et al. (2007); Aguado et al. (2008); Jiang et al. (2008); You et al. (2010). For suitable parameter regimes, the honeycomb model is in the same quantum phase as the toric code in the sense that the two are adiabatically connected, and this guarantees that the ground space properties and quasiparticle excitations of the two models are closely related Kitaev (2006); Hastings and Wen (2005); Hastings (2007); Dusuel et al. (2008); Schmidt et al. (2008); Vidal et al. (2008). In other parameter regimes, the honeycomb model is also a promising avenue for the realization of a non-abelian topological phase, the Ising anyons, which has been proposed as a candidate for topological quantum computing Bravyi (2006); Freedman et al. (2006). Moreover, while the special properties such as the exact degeneracy, harmonic spectrum, and RG fixedness mentioned above are all lost in the honeycomb model, it is still exactly solvable Kitaev (2006); Feng et al. (2007); Chen and Nussinov (2008); Kells et al. (2009a); Pedrocchi et al. (2011); Schmoll and Orús (2017). For these reasons it is desirable to understand the performance of the honeycomb model as a topological quantum memory.

Here we study quantum error correction in the honeycomb model using a combination of analytic and numerical methods. Our main results are summarized in the next subsection. When considering the honeycomb model rather than the toric code, the behavior as a quantum memory will be affected by several distinct effects. Firstly, though the ground spaces of the honeycomb model and the toric code where the quantum information is stored are related, they are not identical, which will lead to different responses to local operations or perturbations. Additionally, on finite system sizes, the code space of the honeycomb model is not precisely degenerate Kells et al. (2008, 2009b), which will result in dephasing effects between encoded states. Another significant source of distinction between the two systems is that we expect them to be subject to different noise channels in the lab. Topological quantum memories are typically assumed to be exposed to some kind of thermal noise source, and so the details of the noise processes will depend heavily on the structure of the spectrum of the Hamiltonian.

We expect the phenomena seen in the honeycomb model to be representative of more general topological phases away from the finely tuned models typically studied. While many properties of topological phases are universal, the details and performance of error correction protocols can be expected to vary as we move around a given topological phase. Existing studies of topological error correction Dennis et al. (2002); Duclos-Cianci and Poulin (2010a, b); Wang et al. (2010a, b); Duclos-Cianci and Poulin (2014); Bravyi and Haah (2013); Bombin et al. (2012); Wootton and Loss (2012); Anwar et al. (2014); Watson et al. (2015); Hutter et al. (2014); Bravyi et al. (2014); Wootton (2015); Fowler (2015); Andrist et al. (2015); Brell et al. (2014); Wootton et al. (2014); Duclos-Cianci and Poulin (2013); Hutter et al. (2015); Hutter and Wootton (2016); Wootton and Hutter (2016); Burton et al. (2017); Herold et al. (2017, 2015); Dauphinais and Poulin (2017); Torlai and Melko (2017); Darmawan and Poulin (2017) are almost exclusively focused on particular representatives of a given phase, usually the RG fixed point, where in 2D errors have a natural interpretation as precisely localized quasi-particles that can be analyzed using the tools of topological quantum field theory. Here we analyze some of the relevant phenomenology associated with being at a point in a topologically ordered phase that is not an RG fixed point, and in particular the effect of these phenomena on quantum information storage protocols. The advantage of studying such phenomena within the honeycomb model is that it retains some tractability that will be lost in the general case, and so at least some additional insights can be gained without appealing to challenging simulations of general quantum many-body systems.

Summary of results

Our study of quantum error correction in the honeycomb model divides the novel phenomenology into two broad classes: coherent and incoherent effects. We treat these separately using analytic methods for the coherent errors in Sec. III and numerical methods for the incoherent errors in Sec. IV.

Coherent effects are those related to the approximate degeneracy of the ground space and the approximate local indistinguishability of different code states. We will provide explicit exponential bounds on these approximations. The main theorems can be summarized as follows; see Sec. III for more precise and detailed formulations. Our first result is a bound on the energy splitting among the four nearly degenerate ground states on the honeycomb model with linear size $N$ (i.e. $2N^{2}$ spins total) with toroidal boundary conditions.

Result 1** (Approximate degeneracy, Thm. 9).**

The energy splitting $\Delta E$ between any pair of ground states in the four-fold ground space of the honeycomb model is bounded by

[TABLE]

where $N$ is the linear size of the lattice, and $c_{1},a>0$ are constants determined by the coupling coefficients of the honeycomb Hamiltonian.

Thus, the dephasing rate between superpositions of ground states is exponentially small in the system size. It is also worth noting that the constant $a$ appearing here and elsewhere is strictly positive in the entire gapped phase of the honeycomb model. Our next result bounds the difference in expectation values of operators supported on sufficiently small regions.

Result 2** (Approximate local indistinguishability, Thm. 16).**

For an arbitrary operator $O_{A}$ located within any simply connected local region $A$ of size $|A|\leq c_{2}N/\log N$, the difference in expectation value between any pair of ground states $|{\psi_{1}}\rangle$ and $|{\psi_{2}}\rangle$ of the honeycomb model is bounded as

[TABLE]

for constants $c_{2},c_{3},a>0$ determined by the coupling coefficients.

Intuitively, this suggests that the honeycomb model ground space is a quantum code capable of correcting errors in these regions $A$. This is indeed the case, and our next result allows us to bound the error $\epsilon(N)$ of a recovery map applied to the code space where the error is measured in the Bures distance; see Definition 18.

Result 3** (Approximate correctability, Thm. 20).**

For any noise channel supported in any simply connected region $A$ with size $|A|\leq c_{4}N/\log N$, the code space is $\epsilon(N)$-correctable in the Bures distance with

[TABLE]

for constants $c_{4},c_{5},a>0$ determined by the coupling coefficients.

We provide explicit expressions for the constant $a$ in Eq. (41) and for $c_{1},\ldots,c_{5}$ in the formal statements of these theorems, Thm. 9, Thm. 16, and Thm. 20 respectively.

The second class of phenomena we study arises from incoherent effects related to the interaction of the system with a thermalizing environment. In Sec. IV we employ numerical tools to study these effects and to evaluate the performance of a standard error correction protocol in their presence, assuming ideal syndrome measurements. Our results can be summarized as follows:

• •

In the regimes studied, the appearance of non-topological excitations does not qualitatively affect the error correction properties.

• •

Asymptotically, the memory lifetime scales exponentially in inverse temperature, as expected Brown et al. (2016).

• •

In the low-temperature finite-size regime, where the diffusion of a single pair of topological excitations is the dominant failure mode, we anticipate an improved scaling of the memory lifetime with system size as compared to the toric code.

This final effect is a result of the structure of the honeycomb Hamiltonian energetically suppressing some local errors more than others. In the low-temperature, finite-size, perturbative regime we probe, with strictly local couplings to the thermal bath, this leads the topological excitations to be constrained to 1D motion instead of the 2D motion of the standard toric code Freeman et al. (2014); Brown et al. (2016). Though this effect relies on the special structure of this regime, it may be exploitable in suitable experimental implementations of a topological quantum memory.

The remainder of the paper is organized as follows. We review the honeycomb model and establish notation, terminology, and some known results in Sec. II. Then Sec. III reviews known results from the general theory of topological phases before proving the three Results listed above. In Sec. IV we discuss our numerical simulation model for quantum error correction and the results of our simulations. We conclude in Sec. V.

II The honeycomb model

The Kitaev honeycomb model is defined with qubits residing on the vertices of a honeycomb lattice, as shown in Fig. 1. In our treatment, the vertical edges will be preferred, and it will be convenient to consider the (square) lattice $\Lambda$ formed by placing a site at each vertical edge. These vertical edges will often be referred to as “dimers”, motivated by the perturbative treatment of this model sketched below Dusuel et al. (2008); Schmidt et al. (2008); Vidal et al. (2008). At each site $\mathbf{q}=(q_{x},q_{y})$ of $\Lambda$, there are two qubits, corresponding to the lower and upper ends of the dimer (denoted $\bullet$ and $\circ$ respectively). We also choose the boundary conditions of $\Lambda$ to be periodic, so that the edges are identified as shown in Fig. 1, and denote its linear dimensions by $N$ (taken to be even for convenience). The Hamiltonian of the model is given by

[TABLE]

for some coupling constants $J_{\alpha}>0\;$ ($\alpha\in\{x,y,z\}$), and with

[TABLE]

where $\sigma^{\alpha}_{i}$ is the relevant Pauli operator at site $i$. We will call the $K^{\alpha}_{\mathbf{q}}$ link generators to distinguish them from the more general class of link operators that they generate under multiplication.

II.1 Symmetries

Though the terms in this Hamiltonian do not commute with one another, constants of motion can be constructed as products of these terms. In particular, for each site $\mathbf{q}$ of $\Lambda$ we can define an operator

[TABLE]

such that $[H,W_{\mathbf{q}}]=0$ for all $\mathbf{q}$ and $[W_{\mathbf{q}},W_{\mathbf{q^{\prime}}}]=0$ for all $\mathbf{q},\mathbf{q^{\prime}}$. These $W_{\mathbf{q}}$ are naturally associated to the hexagonal plaquette of the honeycomb lattice above site $\mathbf{q}$. Two other privileged operators can also be defined as

[TABLE]

Again, it can easily be verified that these operators commute with $H$, with each other, and with all $W_{\mathbf{q}}$. The honeycomb model can be considered as a subsystem code with gauge generators identified with terms in the Hamiltonian (i.e. the link generators), and in this setting the stabilizer group of the code is generated by the $W_{\mathbf{q}}$, $L_{x}$, and $L_{y}$, while the logical algebra is trivial Suchara et al. (2011). However, this is not the interpretation we will use of this model as a quantum code.

The honeycomb model can be in several phases for different regimes of the parameters $J_{x},J_{y},J_{z}$. In this work, we will consider the regime

[TABLE]

where the model is in the same phase as the toric code Kitaev (2006). We will often refer to this as the “abelian” or “gapped” phase in order to distinguish it from the qualitatively different gapless non-abelian phase that can be found in a different parameter regime of the model. For choices of parameters that lead to a gap, the conditions (12) can be imposed without loss of generality. In the abelian phase, there is a 4-dimensional ground space that is separated from other Hamiltonian eigenspaces by a constant gap independent of the system size. As $N$ diverges, this ground space will converge to an exact degeneracy of its energy eigenvalues. A basis for this space can be labeled by $\pm 1$ eigenvalues of the $L_{x}$ and $L_{y}$ operators, and all such states will be $+1$ eigenstates of all $W_{\mathbf{q}}$ operators. When considering the honeycomb model in this phase as a quantum memory, we will refer to the four-fold ground space as the code space $\mathcal{C}_{\mathrm{khm}}$, and as such the $L_{x}$ and $L_{y}$ operators will be logical operators for this code. Eigenstates of $L_{x}$ and $L_{y}$ within the code space will be called code basis states, and be labeled by $\ell=(l_{x},l_{y})$ for $l_{x},l_{y}=\pm 1$ the eigenvalues of $L_{x}$ and $L_{y}$ respectively.

II.2 Exact ground states

While the honeycomb model was originally solved by Majorana fermionization Kitaev (2006), it can also be solved in several other ways Feng et al. (2007); Chen and Nussinov (2008); Kells et al. (2009a); Pedrocchi et al. (2011); Schmoll and Orús (2017). The essential ingredient for our analysis of coherent noise follows the treatment of Ref. Kells et al. (2009a), which solves the honeycomb model without using either Majorana fermions or redundant degrees of freedom. The starting point of this analysis is to describe the model in terms of hard-core bosons and effective spins, followed by fermionization of the hard-core boson degrees of freedom. Following common terminology, we will refer to $-1$ eigenstates of a $W_{\mathbf{q}}$ operator as having a vortex at site $\mathbf{q}$, and $-1$ eigenstates of a $K^{z}_{\mathbf{q}}$ operator as having a broken dimer at link $\mathbf{q}$. By projecting into common eigenspaces of all $W_{\mathbf{q}}$, $L_{x}$, and $L_{y}$, we can solve the model sector by sector. Associating Dirac fermion creation operators $c_{\mathbf{q}}^{\dagger}$ to each site of $\Lambda$, we can represent the link generators in the vortex-free sector as

[TABLE]

Within this sector, a Fourier transform yields a quadratic effective Hamiltonian

[TABLE]

with $\mathbf{k}=(k_{x},k_{y})$, and where

[TABLE]

$H_{\mathrm{vf}}$ can then be diagonalized as

[TABLE]

by the Bogoliubov-Valatin transformation

[TABLE]

where

[TABLE]

The ground states are fermionic vacuum states annihilated by $\gamma_{\mathbf{k}}$ for all $\mathbf{k}$, and have energy

[TABLE]

where $\ell$ denotes the corresponding eigenvalues of $L_{x}$ and $L_{y}$, and where

[TABLE]

Since the code basis state for each $\ell$ is the ground state of a free-fermion Hamiltonian, we can easily calculate expectation values in these states.

Under periodic boundary conditions, the energy dependence on $\ell$ can be considered to arise from the periodicity (anti-periodicity) of $c_{\mathbf{q}}$ in Eq. (15) for $l_{\alpha}=-1$ (+1 respectively). This can conveniently be absorbed into the allowed momenta to give Kells et al. (2009a)

[TABLE]

for $\alpha\in\{x,y\}$ and $0\leq n_{\alpha}<N$.

Note that in the thermodynamic limit, $N\rightarrow\infty$,

[TABLE]

and similarly for $k_{y}$. Since the $l_{\alpha}$-dependent term in Eq. (24) vanishes as $N\rightarrow\infty$, the momenta $x$ and $y$ become independent of $\ell$ for all ground states. Thus the honeycomb model has an exact four-fold degeneracy in the thermodynamic limit. In this limit the energy density per site of $\Lambda$ (i.e. per dimer) becomes

[TABLE]

II.3 In the perturbative regime

Though the honeycomb model can be exactly solved, for example as outlined above, in parts of our study it will be more convenient to use a perturbative approach to analyzing the model, following Ref. Vidal et al. (2008). In the regime where $J_{z}\gg J_{x},J_{y}$, the honeycomb model can be considered as a perturbation away from the toric code. In this limit, the unperturbed Hamiltonian is taken as $H_{0}=-J_{z}\sum_{\mathbf{q}}K^{z}_{\mathbf{q}}$ and the perturbation

[TABLE]

is such that $H=H_{0}+V$.

Our analysis of this model makes use of the perturbative continuous unitary transformation (PCUT), which defines an effective Hamiltonian related to the original Hamiltonian by

[TABLE]

such that $H_{\mathrm{eff}}$ commutes with $H_{0}$, and $U\to I$ as $\frac{J_{x}}{J_{z}}$ and $\frac{J_{y}}{J_{z}}$ vanish. In the thermodynamic limit, this can be calculated as Vidal et al. (2008)

[TABLE]

where $E_{0}$, $\mu$, $C$ and $D$ are constants that can be calculated perturbatively, and $R$ is a string operator that commutes with $K^{z}_{\mathbf{q}}$ at each site except $\mathbf{q}^{\prime}_{1}$ and $\mathbf{q}^{\prime}_{n}$. The precise form for the constants and $R$ will not be crucial in our analysis; the key point is that each of the terms in Eq. (29) is a Pauli operator. The interested reader can consult Ref. Vidal et al. (2008) for details.

Notice that Eq. (29) commutes with $L_{x}$ and $L_{y}$, and so is exactly degenerate. For finite size systems, there are also corrections to Eq. (29) that include loop operators $L_{x}$ and $L_{y}$ Kells et al. (2008, 2009b). However, these appear only at order $N_{x}$ and $N_{y}$ respectively in $\frac{J_{x}}{J_{z}}$ and $\frac{J_{y}}{J_{z}}$. The effect of these corrections is to lift the degeneracy between different eigenstates of these operators. As we will show in Sec. III.2, these effects are exponentially small in the system size, and so we will typically neglect them in our treatment of $H_{\mathrm{eff}}$.

The Hamiltonian $H_{\mathrm{eff}}$ is effective in the sense that observables acting on $H_{\mathrm{eff}}$ are related to those on the original system by the transformation $U$. This can also be computed perturbatively if desired, taking local operators to quasi-local operators. For the most part when considering the honeycomb model in the perturbative regime, we will work in the effective picture for simplicity.

When restricted to the lowest energy eigenspace of $H_{0}$ (the no-broken-dimer sector) with corresponding projector $\Pi_{0}$, the effective Hamiltonian simplifies to

[TABLE]

This can be further simplified by considering only the leading non-trivial order terms. These terms appear at fourth order, and yield the Hamiltonian

[TABLE]

Within the image of $\Pi_{0}$ (that is, the common $+1$ eigenspace of all $K^{z}_{\mathbf{q}}$), we can consider a 2-dimensional effective spin at each site $\mathbf{q}$ of the lattice with a computational basis given by eigenstates of the $Z_{\mathbf{q},\bullet}$ operator. Using this identification, the Hamiltonian $\Pi_{0}H_{\mathrm{eff}}^{(4)}\Pi_{0}$ gives the toric code, with the $W_{\mathbf{q}}$ playing the role of the stabilizer operators as was first noted by Kitaev Kitaev (2006).

For this reason, it will often be convenient to characterize excitations of the honeycomb model in the perturbative regime by eigenvalues of the $K^{z}_{\mathbf{q}}$ and $W_{\mathbf{q}}$ operators. Although a common eigenstate of these operators is not generally an eigenstate of either $H$ or $H_{\mathrm{eff}}$, the ground space of $H_{\mathrm{eff}}$ can be simply found as the $+1$ eigenspace of all $K^{z}_{\mathbf{q}}$ and $W_{\mathbf{q}}$ operators. It is thus convenient to consider the abelian stabilizer group generated by the $K^{z}_{\mathbf{q}}$ and $W_{\mathbf{q}}$, and consider $-1$ eigenvalues for any of these operators as localized quasiparticle excitations—broken dimers or vortices respectively. In the perturbative limit, vortices can be considered as toric code quasiparticles, while broken dimers are simply localized non-topological quasiparticles.

III Coherent noise

III.1 Properties of topological phases

Systems within the same topological phase share common properties that are relevant for error correction. In the following subsections, we will make use of the particular algebraic structure of the honeycomb model to extract bounds on its behavior, but it is also possible to derive results based on general topological properties of the system. This will give some idea of the asymptotic behavior of the model, without giving at all tight bounds. In this subsection we will derive several bounds which motivate the more precise study succeeding it.

The strategy in this analysis will be to relate the properties of the honeycomb model to the known properties of the standard toric code, simply through their being in the same topological phase. A topological phase is defined as an equivalence class of (quasi-)local Hamiltonians connected by smooth, uniformly gapped paths Chen et al. (2010). In order to make this concrete, let us explicitly specify a path of Hamiltonians relating the honeycomb code space to the toric code ground space:

[TABLE]

At $H(0)=-\sum_{\mathbf{q}}\left(J_{w}W_{\mathbf{q}}+J_{z}K^{z}_{\mathbf{q}}\right)$ it is clear that we have an encoded version of the toric code. For our choice of boundary conditions and encoding, this has distance $d=N$ Bombin and Martin-Delgado (2007). That is, it is a $[[2N^{2},2,N]]$ code. It satisfies exact local indistinguishability and it has an exactly degenerate code space Kitaev (2003, 2006). At $H(1)=-\sum_{\mathbf{q}}\left(J_{w}W_{\mathbf{q}}+J_{x}K^{x}_{\mathbf{q}}+J_{y}K^{y}_{\mathbf{q}}+J_{z}K^{z}_{\mathbf{q}}\right)$, we have the honeycomb model Hamiltonian plus some additional $J_{w}W_{\mathbf{q}}$ terms. Since these terms commute with the honeycomb Hamiltonian, and it is known that the ground space of the model is in the common $+1$ eigenspace of all $W_{\mathbf{q}}$, it is easy to see that $H(1)$ has the same code space as the honeycomb Hamiltonian for the same $J_{x},J_{y},J_{z}$.

Within the vortex-free sector, the honeycomb model is known to have a gap between the code space and the next lowest eigenvalue given by $\delta E=2(J_{z}-J_{x}-J_{y})$ in the thermodynamic limit Kitaev (2006). Since $H(s)$ has energy at least $4J_{w}$ for any state orthogonal to the vortex-free sector, by choosing $J_{w}=\frac{1}{4}\delta E$, we can guarantee that $H(s)$ also has a uniform gap above the code space of at least $\delta E$.

The key technical tool underlying the kinds of general topological arguments we will discuss here is the quasi-adiabatic continuation mapping between the code spaces of $H(0)$ and $H(1)$. There are several formulations of the quasi-adiabatic continuation Hastings (2004); Hastings and Wen (2005); Osborne (2007); Nachtergaele and Sims (2007) that would yield slightly different bounds, but for concreteness we follow Ref. Bachmann et al. (2012). In order to be explicit about the locality properties of the quasi-adiabatic continuation, we note that all the Hamiltonians $H(s)$ satisfy a Lieb-Robinson bound Lieb and Robinson (1972), a version of which we state here for completeness.

Lemma 4** (Lieb-Robinson bound Lieb and Robinson (1972)).**

Under time evolution for time $t$, there exist constants $C,\alpha,v>0$ such that observables $O_{A}$ and $O_{B}$ whose support lies on disjoint regions $A,B\subset\Lambda$ satisfy

[TABLE]

for all $s$, where $d(A,B)$ is the distance between regions $A$ and $B$, and $\|\cdot\|$ is the operator norm.

The quasi-adiabatic continuation $U_{\mathrm{qac}}(s)$ maps the code space of $H(0)$ to that of $H(s)$. It is shown in Ref. Bachmann et al. (2012) that a Lieb-Robinson bound implies that the effect of the quasi-adiabatic continuation on an operator whose support is contained within a region $A$ can be approximated well by an operator $V_{A^{r}}(s)$ localized within an $r$-ball around $A$ (we will suppress $s$ labels for clarity).

Lemma 5** (Locality of quasi-adiabatic continuation Hastings (2004), following (Bachmann et al., 2012, Thm. 3.4)).**

The action of the quasi-adiabatic continuation on an operator $O_{A}$ supported on $A$ can be approximated by $V_{A^{r}}$ to an accuracy

[TABLE]

for some constants $\kappa,C^{\prime}$ (depending on coupling constants of the Hamiltonian), and for $\kappa r>c$ for some absolute constant $c$.

Given these lemmas, the next theorem provides a convenient way to summarize their implications for error correction. Following Ref. Flammia et al. (2017), define an $[[n,k,d,\delta,r]]$ code to consist of $k$ logical qubits encoded in $n$ physical qubits, where every region of size less than $d$ is $(\delta,r)$-correctable. A region $A$ is $(\delta,r)$-correctable if an erasure of qubits within $A$ can be corrected to error $\delta$ in the Bures distance, by acting on an $r$-ball enclosing $A$.

Theorem 6** (Smoothness of code parameters (Flammia et al., 2017, Thm. 16)).**

Given a code space $\mathcal{C}_{0}$ with parameters $[[n,k,d]]$ for $d=\Omega(n^{\gamma})$ for some $\gamma>0$, and a unitary $U$ whose action on local regions can be approximated by a local operator as in Eq. (34), $\mathcal{C}_{1}\equiv U\mathcal{C}_{0}U^{\dagger}$ is an $[[n,k,d-2r,\delta(r),r]]$ approximate code, with $\delta(r)\leq 2C^{\prime}\left(\kappa r\right)^{10}\mathrm{e}^{-\frac{2}{7}\frac{\kappa r}{\ln^{2}\kappa r}}$.

According to Thm. 6, we immediately see that the honeycomb model is a $[[2N^{2},2,N-2r,\delta(r),r]]$ approximate error correcting code, for any $r>\frac{c}{\kappa}$. We expect that this bound could be modified in several ways by using alternative formulations of the quasi-adiabatic continuation. In particular, a similar bound should hold for $r\sim 1$.

An intuitive corollary of this theorem is that the honeycomb model is approximately locally indistinguishable. The standard toric code (like other topologically ordered models at the RG fixed point) is well-known to satisfy local indistinguishability in the sense that all codestates look identical when restricted to any local region. In contrast, the states of the honeycomb code space will not exactly satisfy this condition. Thus instead of considering exact local indistinguishability, we consider its approximate version, which we define as follows:

Definition 7** (Approximate local indistinguishability).**

Consider regions $A$ with linear size less than $N^{*}$, and let $O_{\!{A}}$ be any Hermitian operator whose support is contained in $A$, i.e. the linear size of $\mbox{supp}(O_{\!{A}})$ is smaller than $N^{*}$. A subspace $\mathcal{C}$ is $\Delta$-approximately indistinguishable on length scale $N^{*}$ if, for all $A$ and all $O_{\!{A}}$,

[TABLE]

A [math]-approximately indistinguishable space is said to be exactly locally indistinguishable.

A family of subspaces $\{\mathcal{C}_{N}\}_{N}$ of models of lattice size $N$ is said to be approximately locally indistinguishable if there exists $\beta>0$, $N_{0}<\infty$, and superpolynomially decaying function $\Delta$ such that for all $N>N_{0}$, $\mathcal{C}_{N}$ is $\Delta(N)$-approximately indistinguishable on length scale $N^{*}\propto N^{\beta}$. An approximately locally indistinguishable family of states with $\Delta(N)=0$ is said to be exactly locally indistinguishable.

In discussing local indistinguishability, we will typically neglect the distinction between the code space of the honeycomb model on a fixed lattice, and the family of code spaces of the honeycomb model for all lattice sizes. Thus we will abuse the language by saying that the honeycomb code space is approximately locally indistinguishable, it being implicit that the more precise statement is that the family of code spaces of the honeycomb model for all lattice sizes is approximately locally indistinguishable.

Our definition of local indistinguishability differs slightly from related definitions in the literature, since it is typical to either consider only exact local indistinguishability Bravyi et al. (2010); Bravyi and Hastings (2011), or to consider a more restrictive definition of local topological quantum order that includes both local indistinguishability and local consistency Michalakis and Zwolak (2013); Cirac et al. (2013).

In order to show that the honeycomb code space is approximately locally indistinguishable, it will be convenient to recall Thm. 7 of Ref. Flammia et al. (2017), which states that

Theorem 8** (A correctable region avoids logical operators (Flammia et al., 2017, Thm. 7)).**

If region $A$ is $(\delta,r)$-correctable, then for any logical unitary $U_{L}$, there exists an operator $V_{\bar{A}}$ supported only on the complement of $A$, such that

[TABLE]

for $\Pi$ the projector to the code space.

This implies that such a $(\delta,r)$-correctable region $A$ is also $8\sqrt{\delta}$ indistinguishable, since there exists logical unitary $U_{L}$ that maps between any given pair of ground states $U_{L}|\psi_{1}\rangle=|\psi_{2}\rangle$ such that

[TABLE]

where the last line follows from using the above theorem, the triangle inequality, and $[O_{A},V_{\bar{A}}]=0$.

For the honeycomb model, regions $A$ of linear size $N^{*}$ contain at most $2(N^{*})^{2}$ qubits. Thus, for $r<\frac{N}{2}-(N^{*})^{2}$, all such regions satisfy $|A|<d=N-2r$ and are therefore $(\delta(r),r)$-correctable, with $\delta(r)$ as in Thm. 6. By setting $N^{*}=\frac{1}{2}N^{\frac{1}{2}}$ and $r=\frac{N}{2}-(N^{*})^{2}-1$, we see that $8\sqrt{\delta(r)}$ is superpolynomially decaying in $N$. Note that this only holds when $r$ is large enough that Lemma 5 applies, but this will be true for all $N>N_{0}=\frac{4c}{\kappa}+4$. Then Eq. (40) immediately implies that the honeycomb code space is approximately locally indistinguishable.

Alternatively, an analogous result could have been derived directly by considering local operators of the toric code, for which exact local indistinguishability holds, evolved under the truncated quasi-adiabatic continuation $V_{A^{r}}(s)$, and using the bound Eq. (34).

A further implication is that the splitting of the code space of the honeycomb model is superpolynomially suppressed in system size. This can be seen by considering the expectation value of the Hamiltonian in an arbitrary ground state, noting that any local operator can only distinguish the code states by a superpolynomially small amount in system size, and the Hamiltonian is composed of a polynomial number of local operators. This gives that the energies of different code states can only differ by an superpolynomially small amount.

These general properties of approximate degeneracy, approximate local indistinguishability, and approximate code distance motivate the study of the honeycomb model for quantum information storage. By making use of specific properties of the honeycomb model, as opposed to general properties of the toric code phase, in the following sections we will be able to rigorously make much more precise statements about the behavior of this system.

III.2 Approximate degeneracy

In this section, we consider more carefully the approximate degeneracy of the honeycomb Hamiltonian. In the context of error correction, this corresponds to the resilience of information stored in the honeycomb code space against dephasing errors. An intrinsic dephasing effect exists in finite size systems even without the presence of other external noise, due to the inexact degeneracy of the code space.

This approximate degeneracy has been previously studied in Ref. Kells et al. (2008, 2009b), though their perturbative analysis applies only in the limit that $J_{z}\gg J_{x},J_{y}$. In the previous section, we were only able to argue on general grounds that the splitting asymptotically decays superpolynomially. In contrast, here we rigorously and precisely quantify the splitting of ground state energy within the whole abelian phase.

Recall that without loss of generality in the abelian phase, we take $J_{z}>J_{x}\geq J_{y}>0$ and $J_{z}>J_{x}+J_{y}$ as in Eq. (12). Here and elsewhere it will be convenient to introduce a constant

[TABLE]

which is strictly positive everywhere in the gapped phase. The main result of this section is then the following:

Theorem 9**.**

(Energy splitting bound) Let $\Delta E=\left|{E_{\ell}-E_{\ell^{\prime}}}\right|$ be the energy splitting between any pair of (pseudo-ground) states in the four-dimensional code space of the honeycomb model with coupling constants as in Eq. (12). Then

[TABLE]

where $a>0$ is given by a=\ln\bigl{(}\tfrac{J_{z}-J_{y}}{J_{x}}\bigr{)}.

This result shows that, in a superposition of code states, dephasing due to energy splitting decreases exponentially in the linear system size within the entire gapped phase. Moreover, this tendency vanishes exactly where the gap closes on the phase boundary $J_{z}=J_{x}+J_{y}$, where $a\rightarrow 0$.

Proof. Our strategy will be to look at the thermodynamic limit of the model and show that finite-size approximations to this limit converge to the asymptotic result quickly in the linear size of the system. Rather than bound the finite-size corrections perturbatively, we will use arguments from approximation theory that give exponential accuracy.

The key ingredient of our proof is a bound on the error of the rectangular rule from numerical integration. Consider a function $f:\mathbb{R}\rightarrow\mathbb{R}$, and define the error for the $N$-point rectangular rule by

[TABLE]

Suppose that $f$ is also analytic and $2\pi$-periodic. Then there exists a strip in the complex plane $D_{0}=\mathbb{R}\times(-a_{0},a_{0})\in\mathbb{C}$ with $a_{0}>0$ such that $f$ can be extended to a holomorphic and $2\pi$-periodic bounded function $f:D_{0}\rightarrow\mathbb{C}$. We will choose the largest possible strip, so that $a_{0}$ is defined as the supremum of the half-width over of all such strips for $f$. Let $Q_{0}$ be an upper bound for $\left|{f}\right|$ on $D_{0}$; in particular we choose $Q_{0}=\sup_{D_{0}}\left|{f}\right|$. Then by a theorem of Davis Davis (1959) the error for the rectangular rule can be bounded by

[TABLE]

See (Kress, 1998, Thm. 9.28) for a simple proof.

For our purposes, the function $f$ of interest will be the energy density of any of the four honeycomb code basis states. In order to use Davis’ theorem, we must first extend it to apply to a two-dimensional function $f(x,y)$. We define the error of two-dimensional integration

[TABLE]

which can be expressed in terms of the errors of one-dimensional integrations in $x$ and $y$:

[TABLE]

If for all real values $y$, $f(x,y)$ is analytic and $2\pi$-periodic in $x$, then for every value of $y$ there exists a strip of width $a(y)>0$, D_{a}=\mathbb{R}\times\bigl{(}-a(y),a(y)\bigr{)}\in\mathbb{C} and an upper bound $Q(y)=\sup_{x\in D_{a}}|f(x,y)|$ such that

[TABLE]

We also introduce analogous constants $b(x)$ on a strip $D_{b}$ and $R(x)=\sup_{y\in D_{b}}|f(x,y)|$ for the case of the $y$ variable, and we obtain

[TABLE]

According to the definition of $\epsilon_{1}(N;y)$ and $\delta_{1}(N;x)$ we know that

[TABLE]

Using the triangle inequality then gives us an error bound for the two-dimensional rectangular rule,

[TABLE]

Here the optimization is over real values of the periodic variables $x$ or $y$, so a maximum (as opposed to supremum) is appropriate due to compactness. Using the bound from the one-dimensional case, we obtain

[TABLE]

where in the last line we have defined

[TABLE]

The “$\min$” in the definitions of $Q$ and $R$ is not a typo, and the stronger inequality using a minimum instead of a naive maximum follows from the results of Lemma 10 below.

Now let us consider the energy density per dimer of the system in the code basis state $\ell$ of Eq. (22), which can be considered as an approximate evaluation of Eq. (26) as

[TABLE]

where $\epsilon_{\ell}(N)$ is the error term for approximating the two-dimensional integral by a sum for ground state $\ell$. Note that the continuum integral is independent of $\ell$. Then from the triangle inequality we have

[TABLE]

Recall that the mode energy $E_{\mathbf{k}}$ in Eq. (23) is composed of trigonometric functions, is $2\pi$-periodic and is analytic in the gapped phase, so the mode energy satisfies the conditions of Davis’ theorem for both $\ell$ and $\ell^{\prime}$. Since $\epsilon_{\ell}(N)$ and $\epsilon_{\ell^{\prime}}(N)$ are the errors for the same integral, they share the same bound, and

[TABLE]

which follows immediately from the two-dimensional error bound we derived above.

For the coupling constants as in Eq. (12), we will prove in Lemma 10 below that $Q+R\leq 4\sqrt{2}J_{z}$ and b\geq a\geq\ln\bigl{(}\frac{J_{z}-J_{y}}{J_{x}}\bigr{)}. The theorem then follows immediately.  $\Box$

The following lemma provides the explicit values of the constants used in the proof of Thm. 9.

Lemma 10**.**

(Bounds on constants $a,b,Q,R$) For the honeycomb model with coupling constants as in Eq. (12), for all $y$ let $a(y)$ be the supremum half-width over all strips $D_{a}=\mathbb{R}\times(-a(y),a(y))\in\mathbb{C}$ containing an analytic extension in $x$ of $E(x,y)$, let $a=\min_{y}a(y)$, and let $Q=\min_{y}Q(y)=\min_{y}\sup_{D_{a}}\left|{E(x,y)}\right|$. Define $b$ and $R$ in analogy with $a$ and $Q$, but with $x\leftrightarrow y$. Then we have the following inequalities

[TABLE]

as well as

[TABLE]

Proof. To compute the explicit bounds, we have to extend the variables $x$ and $y$ into the complex plane separately. We will focus on the case that $x=x_{\!{R}}+ix_{\!{I}}$ and $y\in\mathbb{R}$, and since $x$ and $y$ are symmetric in the function we consider, the bounds in two directions have the same form with $x$ and $y$ exchanged. At all times we restrict our discussion to the gapped phase with coupling constants as in Eq. (12).

The exponent $a$ in the error bound is determined by the width of the strip where the energy function $E=E(x,y)$ is analytic in the complex plane, so the boundary of the region in $\mathbb{C}$ where $E$ is analytic will give us the width of this strip. Since the square root of a complex number is analytic except at its branch cut, the non-analytic region of $E$ satisfies

[TABLE]

By doing the complex extension explicitly with the energy function from Eq. (23), we have the following expressions

[TABLE]

where $A=J_{y}\sin(y)$, $B=J_{y}\cos(y)+　J_{z}$ and $C=\sqrt{A^{2}+　B^{2}}$ . This expression simplifies if we define unit vectors

[TABLE]

such that

[TABLE]

From Eq. (59) we find that for non-analytic $E$, either $\sinh(x_{\!{I}})=0$ or $\hat{c}\cdot\hat{z}^{\perp}=0$, since $C>0$ and $J_{x}>0$. The first case gives $x_{\!{I}}=0$ which is trivial, and the second reveals that $\hat{c}\cdot\hat{z}=\pm 1$. Hence Eq. (59) gives

[TABLE]

There is no solution for $x_{\!{I}}$ in the first case since $-\frac{J_{x}^{2}+C^{2}}{2J_{x}C}<0$. In the second case, we have

[TABLE]

and this inequality reveals the width of strip, $a(y)$. We can make one step further since $\cosh^{-1}(u)=\ln(u+\sqrt{u^{2}-1})$ for $u\geq 1$. The condition is always satisfied because $J_{x},C>0$ in the gapped phase. hence the width of the strip in the gapped phase (where $C^{2}>J_{x}^{2}$) becomes

[TABLE]

and minimizing over the real variable $y$ we find

[TABLE]

The value for $b$ follows the same argument, and yields b=\ln\bigl{(}\frac{J_{z}-J_{x}}{J_{y}}\bigr{)}. The inequality $b\geq a$ follows from elementary algebra using the conditions in Eq. (12).

To compute $Q$, we apply the maximal modulus principle from complex analysis that states that the maximum modulus of a function $E$ which is analytic in an open subset $D$ of the complex plane lies on the boundary $\partial D$ of $\bar{D}$, the closure of $D$. To make the calculation simpler, we take the fourth power of the energy function to eliminate the square root, such that

[TABLE]

The boundary of the strip is along the line x_{\!{I}}=\cosh^{-1}\biggl{(}\frac{J_{x}^{2}+C^{2}}{2J_{x}C}\biggr{)}, and on the boundary the energy function takes the values

[TABLE]

Since only $\hat{z}$ and $\hat{z}^{\perp}$ depend on $x_{\!{R}}$, and they only appear in the term inside an inner product with $\hat{c}$, we can replace the $x_{\!{R}}$-dependence by an angle $\theta$ such that $\hat{c}\cdot\hat{z}=\cos(\theta)$ and $\hat{c}\cdot\hat{z}^{\perp}=\pm\sin(\theta)$. Now we arrive at a simple form for the energy function on the boundary,

[TABLE]

The extreme values of this function are at points where $\tfrac{\partial}{\partial\theta}\left|{E}\right|^{4}=0$, which is equivalent to

[TABLE]

There is no solution for $\theta$ in the $\cos(\theta)$ equation, since for a gapped system (Eq. (12)) we have \bigl{(}C^{2}+J_{x}^{2}\bigr{)}^{2}-4C^{2}J_{x}^{2}>0. Therefore the maximum of the energy function is given by

[TABLE]

and

[TABLE]

Minimizing this expression over the real variable $y$ yields

[TABLE]

The derivation for the constant $R$ is the same, except exchanging $x$ and $y$.

Having derived the expressions for $Q(y)$ and $a(y)$ (as well as their minimums), we now turn to the claimed upper bound

[TABLE]

First let

[TABLE]

Then the maximization on the left hand side can be rewritten in this new variable as simply

[TABLE]

The stationary points of g\bigl{(}u(y)\bigr{)} are where $\frac{\mathrm{d}g}{\mathrm{d}u}\frac{\mathrm{d}u}{\mathrm{d}y}=0$, and $g(u)$ is a strictly decreasing function of $u$ for $u>1$, so there are no stationary points where $\frac{\mathrm{d}g}{\mathrm{d}u}$ vanishes. Therefore the maximum must occur where $\frac{\mathrm{d}u}{\mathrm{d}y}=0$, which simplifies to just $\sin(y)=0$. Again because $g(u)$ is strictly decreasing the maximum occurs when $u$ is a small as possible, and it follows that $y=\pi$ is always the maximizer, justifying our choice to minimize in the definitions of $Q$ and $a$. This establishes the bound in Eq. (82). The analogous inequality involving $R$ and $b$ follows the same argument.

Finally, the upper bound on $Q+R$ then comes from maximizing the coupling constants $J_{x},J_{y}$ in this sum over the region specified by Eq. (12), with the maximum occurring in the limit that $J_{x},J_{y}\to 0$.  $\Box$

III.3 Approximate local indistinguishability

In addition to the intrinsic dephasing generated by finite size effects, another source of coherent noise for quantum information stored in the honeycomb model will be from local perturbations or local operations, representing corrections to the model or experimental imperfections. To analyze the stability of this system against these sources of error, we study the local indistinguishability of the honeycomb code space Kitaev (2003); Bravyi et al. (2010); Bravyi and Hastings (2011); Michalakis and Zwolak (2013).

As argued in Sec. III.1, the honeycomb code space will not be exactly locally indistinguishable, but will satisfy the approximate local indistinguishability condition given in Definition 7. Although this has already been shown, the decay functions and constants in the argument are far from optimal. In this section, we will reprove the approximate local indistinguishability of the honeycomb model in a more precise form. The explicit theorem will be presented as Thm. 16, but before proving this we will introduce several relevant lemmas.

As in Sec. III.2, the key technical ingredient of the proof is a bound on the discrepancy between quantities on finite size lattices and quantities in the thermodynamic limit. Our argument begins by classifying Pauli operators into two sets, the centralizer of the set of $W_{\mathbf{q}}$ and its complement. The complement is shown to satisfy exact local indistinguishability in Lemma 11, which leads us to focus on the set of local Pauli operators that commute with all $W_{\mathbf{q}}$. This set is proved to be the group of link operators (those operators generated by link generators) in Lemma 12. Analogously to the proof of Theorem 9, we can express the expectation value of these link operators as discrete Fourier transforms, and approximate them by continuous Fourier transforms for large $N$. We bound the error of this approximation for string-like link operators in Lemma 14, and general link operators in Lemma 15. Finally in Thm. 16 we express an arbitrary local operator in the Pauli basis to demonstrate its approximate local indistinguishability.

For the remainder of this section, we will consider only regions $A$ formed by a simply-connected set of honeycomb lattice plaquettes. We define the corresponding set of plaquette operators contained in $A$ as

[TABLE]

It is clear that the honeycomb codespace $\mathcal{C}_{\mathrm{khm}}$ possesses the following property:

Lemma 11**.**

For any Pauli $P$ with $P\notin\mathcal{Z}(\mathcal{W}(\Lambda))$ (where $\mathcal{Z}(\mathcal{S})$ denotes the centralizer of set $\mathcal{S}$, and recall that $\Lambda$ is the entire lattice), and all $|{\psi}\rangle\in\mathcal{C}_{\mathrm{khm}}$,

[TABLE]

Proof. It is known that $\mathcal{C}_{\mathrm{khm}}$ is in the vortex-free sector Lieb (1994), and so any $P\notin\mathcal{Z}(\mathcal{W}(\Lambda))$ will map $\mathcal{C}_{\mathrm{khm}}$ into the orthogonal complement of $\mathcal{C}_{\mathrm{khm}}$. Thus $\left\langle{P}\right\rangle_{\psi}=0$ for any $|{\psi}\rangle\in\mathcal{C}_{\mathrm{khm}}$.  $\Box$

Given this observation, in order to prove indistinguishability we can restrict to considering Paulis $P\in\mathcal{Z}(\mathcal{W}(\Lambda))$ (since the Pauli operators form a basis).

A vertex is said to be in the boundary of $A$ if $A$ does not contain all of its neighbors. A link is said to be across the boundary of $A$ if it contains one vertex within $A$ and another vertex outside of $A$. As well as $\mathcal{W}(A)$, the set of plaquette operators within $A$, we define the set of link generators within $A$, and the restrictions of link generators across the boundary of $A$ to $A$ itself as

[TABLE]

Lemma 12** (Classification of Pauli operators).**

Every Pauli $P_{\!{A}}\in\mathcal{Z}(\mathcal{W}(\Lambda))$ with support contained within region $A$ (formed of simply-connected plaquettes) is a link operator in $A$. In other words,

[TABLE]

where $\langle\mathcal{S}\rangle$ denotes the group generated by set $\mathcal{S}$. Further, $\langle\mathcal{K}(A)\rangle$ contains at most $2^{\frac{3}{2}\left|{A}\right|}$ linearly independent elements.

Note that, to avoid confusion between notation for group generation and expectation value, expectation values will always be labelled with a state subscript as in $\left\langle{O}\right\rangle_{\psi}$.

Proof. The number of vertices, edges, plaquettes, and boundary vertices in $A$ are given by $|A|$, $|\mathcal{K}|$, $|\mathcal{W}|$, and $|\mathcal{B}|$ respectively (in this proof we will suppress $A$ arguments of the sets $\mathcal{K}(A)$, $\mathcal{W}(A)$, and $\mathcal{B}(A)$ for clarity). Since $A$ is simply connected, we have from Euler’s formula

[TABLE]

All boundary vertices of $A$ are contained in two edges of $A$ while all other vertices of $A$ are contained in three edges of $A$. Since each edge contains two vertices, by adding up all the edges at each vertex, we have counted every edge exactly twice, giving $2|\mathcal{B}|+3(|A|-|\mathcal{B}|)=2|\mathcal{K}|$. Rearranging, we find

[TABLE]

In general, given a set $\mathcal{S}$ of multiplicatively independent commuting Pauli stabilizers on $n$ qubits (such that $\left\langle{\mathcal{S}}\right\rangle\not\owns-I$), a minimal generator of the centralizer $\mathcal{Z}(\mathcal{S})$ has size $2n-|\mathcal{S}|$ Nielsen and Chuang (2000). We can use this fact to show that $\bar{\mathcal{K}}\equiv\mathcal{K}\cup\mathcal{B}$ generates $\mathcal{Z}_{A}(\mathcal{W})$ (the centralizer of $\mathcal{W}$ in the Pauli group on qubits in $A$) by a simple counting argument.

There are $|\bar{\mathcal{K}}|=|\mathcal{K}|+|\mathcal{B}|$ operators in this generating set, but only $|\bar{\mathcal{K}}|-1$ of them are (multiplicatively) independent, since $\prod_{P\in\bar{\mathcal{K}}}P=I$. We know that $\mathcal{Z}_{A}(\mathcal{W})$ has a minimal generator of size $2|A|-|\mathcal{W}|=3|A|-|\mathcal{K}|-1$. Since $\bar{\mathcal{K}}$ also contains $|\mathcal{K}|+|\mathcal{B}|-1=3|A|-|\mathcal{K}|-1$ independent generators, each of which lies in $\mathcal{Z}_{A}(\mathcal{W})$, we conclude that $\langle\bar{\mathcal{K}}\rangle=\mathcal{Z}_{A}(\mathcal{W})$ as claimed.

It is also easy to see that each operator in $\mathcal{B}$ will anticommute with at least one $W$ stabilizer not contained in $A$, as will all of the operators generated by $\mathcal{B}$, except for the identity $I$ and $\prod_{P\in\mathcal{B}}P$. Since all elements of $\mathcal{K}$ are in $\mathcal{Z}(\mathcal{W}(\Lambda))$, and $\prod_{P\in\mathcal{B}}P=\prod_{P\in\mathcal{K}}P$ can be generated by the elements of $\mathcal{K}$, we conclude that $\mathcal{Z}_{A}(\mathcal{W}(\Lambda))=\langle\mathcal{K}\rangle$.

Since each element of $\langle\mathcal{K}\rangle$ is self inverse and commutes with every other up to a phase, and since its minimal generator has size $|\mathcal{K}|$, the number of linearly independent elements in $\langle\mathcal{K}\rangle$ is simply $2^{|\mathcal{K}|}$. Noting that each edge involves 2 qubits, and every qubit is involved in either 2 or three edges within $A$, it is easy to see that $|A|\leq|\mathcal{K}|\leq\frac{3}{2}|A|$. This bounds the number of linearly independent elements of $\langle\mathcal{K}\rangle$ to be less than $2^{\frac{3}{2}\left|{A}\right|}$.  $\Box$

The result of Lemma 12 allows us to focus on link operators in $\left\langle{\mathcal{K}(A)}\right\rangle$. In order to analyze these operators, it will be convenient to make the following definitions:

Definition 13**.**

An endpoint of a link operator is a qubit with an odd number of link generators acting on it. A link operator is string-like if it has precisely two endpoints.

In fact, we will now see that the expectation value of a link operator in a fixed honeycomb model code basis state depends (up to a phase) only on its endpoints.

Noting that any link operator has an even number of endpoints, we denote the $2M$ endpoints of a link operator $K\in\left\langle{\mathcal{K}(A)}\right\rangle$ as $\mathbf{Q}_{K}=\left\{(\mathbf{q}_{m},s_{m})\right\}_{m=1}^{2M}$, where $s_{m}\in\{+1,-1\}$ corresponding to $\bullet$ and $\circ$ sites, respectively.

Denote the code basis state $\ell=(l_{x},l_{y})$ (the state in the codespace that is also the $l_{x}$ eigenstate of $L_{x}$ and the $l_{y}$ eigenstate of $L_{y}$) by $|{\phi_{\ell}}\rangle$. In any such state, link generators $K^{\alpha}_{\mathbf{q}}$ can be fermionized as in Eq. (II.2). Since $(c^{\dagger}_{\mathbf{q}}\pm c_{\mathbf{q}})^{2}=\pm 1$ and $K^{z}_{\mathbf{q}}=(c^{\dagger}_{\mathbf{q}}+c_{\mathbf{q}})(c^{\dagger}_{\mathbf{q}}-c_{\mathbf{q}})$, we can then calculate

[TABLE]

where $\sim$ reflects neglecting the possible $\pm 1$, $\pm i$ phases. Note that this decomposition is contingent on the fact that $A$ is simply connected, since we can then without loss of generality assume that $A$ does not cross the periodic boundary of the lattice.

By the Fourier transformation in Eq. (15) and the Bogoliubov transformation in Eq. (19), we know that in code basis states,

[TABLE]

This gives for the expectation value of a link operator $K$:

[TABLE]

where

[TABLE]

Since code basis states are ground states of a free fermion Hamiltonian, we can apply Wick’s theorem Wick (1950), and we find

[TABLE]

where the sum is over all involutive derangements $\mu$ of $\mathbf{Q}_{K}$ (i.e. permutations that give rise to valid pairings, satisfying $\mu=\mu^{-1}$ and $\mu(m)\neq m$). Here $\pi_{\mu}$ is a permutation that takes the modes from normal order into an order in which $m$ and $\mu(m)$ are adjacent, and $m$ precedes (follows) $\mu(m)$ for $m<\mu(m)$ ($m>\mu(m)$). Each of the $\frac{(2M)!}{2^{M}M!}$ terms in this expansion has exactly $M$ factors.

Since code basis states $\phi_{\ell}$ are fermion vacuum states, and so are annihilated by $\gamma_{\mathbf{k}_{m}}$ for all $\mathbf{k}_{m}$, we can compute each expectation value explicitly as

[TABLE]

with discrete Fourier coefficients

[TABLE]

Using the expressions (20) and (21), we find

[TABLE]

This allows us to simplify Eq. (96) slightly, since $\left\langle{\nu_{i}\nu_{j}}\right\rangle_{\phi_{\ell}}$ will vanish for $s_{i}=s_{j}$. There are then at most $M!$ non-vanishing terms in this expression (since this is how many pairings of $M$ elements with $M$ elements there are). We denote the sum over these non-vanishing pairings as $\sum_{\tilde{\mu}}$ to distinguish it from the naive sum in Eq. (96).

The reformulation given by equations (96-101) will allow us to bound the difference of expectation values of link operators between any two codestates. Before considering the general case, we will first demonstrate the result for the simpler class of string-like link operators.

Given the general expression of an expectation value for any $K\in\langle\mathcal{K}(A)\rangle$, we bound the difference of expectation values of code basis states, beginning with the following lemma.

Lemma 14** (String-like link operator indistinguishability).**

Let $K\in\langle\mathcal{K}(A)\rangle$ be a string-like link operator with endpoints $\mathbf{Q}_{K}=\{(\mathbf{q}_{1},s_{1}),(\mathbf{q}_{2},s_{2})\}$. Let $\mathbf{\Delta q}=\mathbf{q}_{1}-\mathbf{q}_{2}$, and $\phi_{\ell}$, $\phi_{\ell^{\prime}}$ be any pair of common eigenstates of $L_{x}$ and $L_{y}$ in $\mathcal{C}_{\mathrm{khm}}$ with coupling constants as in Eq. (12). Then for $\|\mathbf{\Delta q}\|_{\infty}\leq\tfrac{N-1}{2}$

[TABLE]

where $c=\frac{4\sqrt{2}}{\sinh(a/2)}$, a=\ln\bigl{(}\tfrac{J_{z}-J_{y}}{J_{x}}\bigr{)}.

Proof. This proof uses the same idea as the proof of Thm. 9, but instead of bounding the difference between the rectangular rule and an integral, we will bound the difference between the discrete Fourier transform (DFT) from its continuous version. Consider a function $f:\mathbb{R}\rightarrow\mathbb{C}$, and define the error for the $N$-point discrete Fourier transform by

[TABLE]

Let us recall the conditions that we used in Thm. 9: Suppose that $f$ is also analytic and $2\pi$-periodic. Then there exists a strip in the complex plane $D_{0}=\mathbb{R}\times(-a_{0},a_{0})\in\mathbb{C}$ with $a_{0}>0$ such that $f$ can be extended to a holomorphic and $2\pi$-periodic bounded function $f:D_{0}\rightarrow\mathbb{C}$. We will choose the largest possible strip, so that $a_{0}$ is defined as the supremum of the half-width over of all such strips for $f$. Let $Q_{0}$ be an upper bound for $\left|{f}\right|$ on $D_{0}$; in particular we choose $Q_{0}=\sup_{D_{0}}\left|{f}\right|$. Our error term for the DFT can be bounded by a theorem of Epstein Epstein (2005) and we find that

[TABLE]

Since our integration is two-dimensional, we define the error of the two-dimensional DFT as

[TABLE]

Following the same procedure as in the proof of Thm. 9, the error bound of two dimensional error based on the one-dimensional case is

[TABLE]

where $a,b,Q,R$ have the same definitions given previously in the approximate degeneracy proof, but of course we will apply them to the function $\hat{h}_{\mathbf{k}}(s_{i},s_{j})$ instead of the energy density.

Now let us consider the expectation values of a string-like link operator in the code basis states $\ell$ and $\ell^{\prime}$. Since there are only two endpoints of a string-like link operator, Eq. (96) consists of a single term, and Eq. (97) allows us to convert this calculation into the approximate evaluation of a continuous Fourier transform with error

[TABLE]

where $\mathbf{x}=(x,y)$. Then from the triangle inequality we have

[TABLE]

Recall that $\hat{h}_{\mathbf{k}}(s_{i},s_{j})$ in Eq. (101) is $\pm 1$ for $s_{i}=s_{j}$ or a function only of basic trigonometric functions of $\mathbf{k}$ for $s_{i}\neq s_{j}$. In the case $s_{1}=s_{2}$, since $\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{\Delta q}}=0$, the expectation values both vanish for all $\mathbf{\Delta q}$, and their difference is 0. In the other case $s_{1}\neq s_{2}$ the integrand $\hat{h}_{\mathbf{k}}(s_{i},s_{j})$ is $2\pi$-periodic, and it is analytic in the gapped phase. For a $2\pi$-periodic function $f(x)$, its Fourier transformation $\int_{0}^{2\pi}f(x+\delta)e^{i(x+\delta)q}\mathrm{d}x=\int_{0}^{2\pi}f(x)e^{ixq}\mathrm{d}x$ for arbitrary real constant $\delta$, so we can apply Epstein’s theorem to bound the error independent of $\ell$, and we find

[TABLE]

where, similarly to the proof of Thm. 9, we define

[TABLE]

and where the implicit function dependence is on $\hat{h}_{\mathbf{k}}(s_{i},s_{j})$. Note that we have chosen to maximize $\tilde{Q}(y)$ and $\tilde{R}(y)$ now because this specific expression does not obviously enjoy the same monotonicity properties that allowed us to use a minimum in the proof of Thm. 9. The constants $a$ and $b$ turn out to be identical to those in the approximate degeneracy bound, since the point of nonanalyticity is still determined by the energy. A similar calculation as in the proof of Lemma 10 then shows that $\tilde{R},\tilde{Q}\leq\sqrt{2}$. We omit the details of these calculations, and the lemma then follows straightforwardly from Eq. (110) with these values for the constants.  $\Box$

Lemma 15** (General link operator indistinguishability).**

Let $K\in\langle\mathcal{K}(A)\rangle$ be a link operator with endpoints $\mathbf{Q}_{K}=\left\{(\mathbf{q}_{m},s_{m})\right\}_{m=1}^{2M}$, $\mathbf{\Delta q}_{mn}=\mathbf{q}_{m}-\mathbf{q}_{n}$, and $\phi_{\ell}$, $\phi_{\ell^{\prime}}$ be any pair of common eigenstates of $L_{x}$ and $L_{y}$ in $\mathcal{C}_{\mathrm{khm}}$ with coupling constants as in Eq. (12). Then for $\max_{m,n}\|\mathbf{\Delta q}_{mn}\|_{\infty}<\tfrac{N-1}{2}$

[TABLE]

where $c=\tfrac{4\sqrt{2}}{\sinh(a/2)}$ and a=\ln\bigl{(}\tfrac{J_{z}-J_{y}}{J_{x}}\bigr{)}

Proof. Given the general form of $K$ in Eq. (96), we first use the triangle inequality to obtain

[TABLE]

Then we use the following identity of a telescoping sum and again apply the triangle inequality,

[TABLE]

where by convention $\prod_{k=1}^{0}=\prod_{l=M+1}^{M}\equiv 1$. Assuming that $\left|{A_{k}}\right|,\left|{B_{k}}\right|\leq 1$, we can write an even simpler inequality,

[TABLE]

Now recall that $\left\langle{\nu_{i}\nu_{j}}\right\rangle_{\phi_{\ell}}$ is the DFT of $\hat{h}_{\mathbf{k}}$, and \bigl{|}\hat{h}_{\mathbf{k}}\bigr{|}=1 in Eq. (101), so we have \bigl{|}\left\langle{\nu_{i}\nu_{j}}\right\rangle_{\phi_{\ell}}\bigr{|}\leq 1. Moreover, from Lemma 14 we have a uniform bound for the difference of a single pair from the string-like case (given the size restriction $\|\mathbf{\Delta q}\|_{\infty}<\tfrac{N-1}{2}$), namely

[TABLE]

Putting these ingredients together, we find that

[TABLE]

where we have used the fact that $M\cdot M!\leq M^{M}=\mathrm{e}^{M\ln M}$.  $\Box$

Theorem 16** (Arbitrary local operator indistinguishability).**

Let $A$ be a simply connected region with volume $\left|{A}\right|\leq\frac{aN/4}{\ln(aN/4+1)}$. Then with coupling constants as in Eq. (12), the code space $\mathcal{C}_{\mathrm{khm}}$ is $\Delta$-approximately indistinguishable on all such regions $A$, i.e.

[TABLE]

for all operators $O_{\!{A}}$ whose support is contained in $A$, and with decay function given by

[TABLE]

Here the constants can be chosen as $c=\tfrac{4\sqrt{2}}{\sinh(a/2)}$ and a=\ln\bigl{(}\tfrac{J_{z}-J_{y}}{J_{x}}\bigr{)}.

For any fixed values of the coupling constants $J_{\alpha}$, this theorem shows that the honeycomb codespace satisfies $\Delta$-approximate local indistinguishability at any length scale N^{*}\leq\mathcal{O}\bigl{(}\sqrt{N/\ln N}\bigr{)}, corresponding to any $\beta<1/2$ in Definition 7, and so is approximately locally indistinguishable. We remark that the proof below does not depend on the geometry of the regions considered, only on their volume and topology (simply connected).

Proof. A general operator supported on $A$ will consist of linear combinations of Pauli operators that either commute or anticommute with the plaquette operators. Operators that anticommute are always trivial in expectation by Lemma 11, so we can expand any nontrivial part of a general operator in terms of $K_{i}\in\langle\mathcal{K}(A)\rangle$ only by Lemma 12. Write a general operator supported on the $K_{i}$ as $O_{\!{A}}=\sum_{i}r_{i}K_{i}$.

Next note that the operators $L_{x}$ and $L_{y}$ commute with the entire link group $\mathcal{K}(A)$. Therefore any observable supported only on elements of the link group breaks into a direct sum labeled by the eigenvalues of the logical operators. In particular, any off-diagonal element of the form $\langle\phi_{\ell}|O_{\!{A}}|\phi_{\ell^{\prime}}\rangle$ must vanish for code basis states $\phi_{\ell}\not=\phi_{\ell^{\prime}}$ (since they are simultaneous eigenstates of $L_{x}$ and $L_{y}$).

Now consider an arbitrary (potentially mixed) state $\rho$ supported only on $\mathcal{C}_{\mathrm{khm}}$. The expectation value of an operator $O_{\!{A}}$ supported on $\mathcal{K}(A)$ is given by $\left\langle{O_{\!{A}}}\right\rangle_{\rho}=\operatorname{Tr}(O_{\!{A}}\rho)$. From the previous paragraph, we know that $O_{\!{A}}$ can be chosen to be diagonal in the basis of states $\{\phi_{\ell}\}$. Therefore by convexity and the variational characterization of eigenvalues, we have

[TABLE]

This implies that the difference between expectation values of any $\rho_{1}$ and $\rho_{2}$ in $\mathcal{C}_{\mathrm{khm}}$ is bounded by

[TABLE]

It then suffices to prove the bound only for these eigenstates.

Using the expansion into Pauli operators inside $A$, we apply Hölder’s inequality to obtain

[TABLE]

Since the $K_{i}$ are all unitary trace-orthogonal Paulis acting on a $2^{2N^{2}}$-dimensional space, we can bound the 2-norm of the coefficient vector $r$ as follows,

[TABLE]

where $\|O_{\!{A}}\|_{\mathrm{F}}$ is Frobenius norm. Because there are at most $2^{3\left|{A}\right|/2}$ linearly independent $K_{i}$ by Lemma 12, we therefore have the bound on the 1-norm

[TABLE]

From the bound in Lemma 15, and the fact that a region with $\left|{A}\right|$ spins has at most $2M\leq\left|{A}\right|$ string endpoints, we obtain

[TABLE]

We have $\left|{A}\right|\geq 6$ since we only consider simply connected plaquettes and the smallest size region must contain at least one plaquette. Then we have the inequality $2^{3\left|{A}\right|/4}\leq\mathrm{e}^{\tfrac{\left|{A}\right|}{2}\ln\tfrac{\left|{A}\right|}{2}}$ and the bound simplifies,

[TABLE]

To achieve our desired decay function with decay rate $a/4$, we require that the region size $\left|{A}\right|$ satisfies $\left|{A}\right|\ln\left|{A}\right|\leq\frac{aN}{4}$. It is straightforward to show that choosing

[TABLE]

gives the desired result for all $a>0$.  $\Box$

For some applications, a slightly different form of this result is desirable that emphasizes the projection into the code space rather than differences of expectation values. We state here an immediate corollary of Thm. 16 that is more convenient in these cases.

Corollary 17**.**

Under the same conditions as Thm. 16, and with $\Pi$ the projector onto the code space $\mathcal{C}_{\mathrm{khm}}$, we have

[TABLE]

Proof. Since $\Pi O_{\!{A}}\Pi$ is codiagonal with $\Pi$, the proof follows immediately from Thm. 16 and the triangle inequality, noting only that $\operatorname{Tr}\Pi=4$.  $\Box$

III.4 Consequences for error correction

In Sec. III.1, local indistinguishability was derived as a consequence of the code properties of the honeycomb model. These two notions are closely connected, and here we will illustrate this by using the indistinguishability results obtained in the previous section to derive consequences for the error correction properties of the honeycomb code.

However, here it will be convenient to use a slightly different formulation of approximate codes to that used previously. In Sec. III.1, we used a notion of approximate codes due to Flammia et. al. Flammia et al. (2017) suited to the topological arguments of that section. This formalism restricts to codes with some notion of geometric locality, and treats erasure errors preferably. In contrast, we will now consider a slightly more abstract formulation of approximate codes due to Beny & Oreshkov Bény and Oreshkov (2010). In this formalism, approximate correctability is defined as follows:

Definition 18**.**

A code defined by an encoding map $\mathcal{E}$ is called $\epsilon$-correctable under a noise channel $\mathcal{N}$ if there exists a recovery map $\mathcal{R}$ such that $d(\mathcal{RNE},I)\leq\epsilon$, where

[TABLE]

is the Bures distance based on the entanglement fidelity, with $\psi$ a purification of $\rho$ and $f(\rho,\sigma)=\operatorname{Tr}\sqrt{\sqrt{\sigma}\rho\sqrt{\sigma}}$ the fidelity.

Ref. Bény and Oreshkov (2010) also provides a useful characterization of codes that are approximately correctable:

Theorem 19** (Approximate code correctability (Bény and Oreshkov, 2010, Cor. 2)).**

A code space $\mathcal{C}$ is $\epsilon$-correctable under a noise channel with elements $\{F_{i}\}$ iff

[TABLE]

for $\Pi_{\mathcal{C}}$ the projector to $\mathcal{C}$, $T_{ij}$ arbitrary matrices, $g_{ij}$ the components of a density operator, and with $d(\Gamma+\mathcal{T},\Gamma)\leq\epsilon$, where $\Gamma(\rho)=\operatorname{Tr}(\rho)\sum_{ij}g_{ij}|{i}\rangle\!\langle{j}|$ and $\mathcal{T}(\rho)=\sum_{ij}\operatorname{Tr}(\rho T_{ij})|{i}\rangle\!\langle{j}|$.

Using Beny & Oreshkov’s conditions together with the result of Corollary 17 yields the following bound on the correctability of the honeycomb code. Recall that as in the previous section, we restrict our attention to simply connected regions composed of whole plaquettes. This easily extends to arbitrary simply connected regions by expanding the region slightly to complete any partial plaquettes.

Theorem 20** (Local approximate correctability of the honeycomb code).**

The honeycomb code space $\mathcal{C}_{\mathrm{khm}}$ is $\epsilon(N)$-correctable under a noise channel with elements supported in any simply connected region $A$ composed of complete plaquettes with volume $|A|\leq\min\{\frac{aN}{8\ln 4},\frac{aN/4}{\ln(aN/4+1)}\}$, where

[TABLE]

and where $c=\tfrac{4\sqrt{2}}{\sinh(a/2)}$ and a=\ln\bigl{(}\tfrac{J_{z}-J_{y}}{J_{x}}\bigr{)}.

Proof. Consider an orthonormal basis for operators supported in $A$ (of size $4^{|A|}$), and take any two basis elements $F^{A}_{i}$ and $F^{A}_{j}$. Then note that $O_{ij}=\left(F^{A}_{i}\right)^{\dagger}F^{A}_{j}$ is also supported in $A$, and so satisfies Eq. (128). This gives that

[TABLE]

for some $T_{ij}=\Pi T_{ij}\Pi$ with $\|T_{ij}\|\leq\tfrac{3}{4}\Delta(N)\|O_{ij}\|\leq\tfrac{3}{4}\Delta(N)$. Setting $g_{ij}=\tfrac{\operatorname{Tr}(\Pi O_{ij})}{\operatorname{Tr}\Pi}$, it can easily be verified that $\Gamma$ is trace-one and Hermitian as required.

Next we relate the fidelity to the trace distance using the Fuchs-van de Graaf inequality Fuchs and van de Graaf (1999), which says that for any pair of states $\rho$ and $\sigma$, we have

[TABLE]

It follows that it suffices to bound

[TABLE]

where $\psi$ is a purification of $\rho$. For any matrix $M$, we have $\|M\|_{1}\leq\dim(M)\max_{i,j}|M_{i,j}|$, and then using the matrix Hölder inequality it follows that for all $\psi$

[TABLE]

Here the second inequality also uses the upper bound on the purification dimension, and subsequent lines use submultiplicativity and Eq. (128).

Putting this together with Eq. (119), we find that

[TABLE]

Since $|A|\leq\frac{aN}{8\ln(4)}$ by assumption, the result follows directly.  $\Box$

This demonstrates that, as with the approximate degeneracy and approximate local indistinguishability of the honeycomb model, its approximate code properties can also be demonstrated more precisely than was shown using general topological arguments of Sec. III.1. However, the present result is slightly weaker in the sense that it only holds for simply-connected regions $A$. We anticipate that this requirement could be relaxed by modifying the arguments in Sec. III.3 if desired, but that this would be relatively involved.

IV Incoherent noise

IV.1 Simulation of thermalization processes

A general framework for simulating many-body quantum systems interacting with a thermalizing environment is to model the dynamics by a continuous-time Markov process, where we assume that each spin couples to an independent environment (for a general discussion on thermalization in these kinds of systems, see e.g. Ref. Brown et al. (2016)). A common form for such dynamics is given by the spin-boson model, which considers errors $F$ that decrease the energy of the system by $\Delta_{F}$. It assigns rates to such processes as

[TABLE]

for inverse temperature $\beta$. These dynamics have desirable properties such as being motivated by a microscopic model, having the Gibbs state as a fixed state, and satisfying the detailed balance condition. To simulate thermalization using this model, a set of (typically ergodic) allowed errors is chosen, and processes are chosen at random according to their relative rates for a Poisson distributed number of time steps around the desired simulation time. For more details about these standard simulation techniques, see Chesi et al. (2010); Brown et al. (2016).

In order for this approach to be well motivated, it is important that an error $F$ will map energy eigenstates to energy eigenstates, so that $\Delta_{F}$ is well-defined. This is the case for Pauli errors $F$ and stabilizer Hamiltonians, but in general, Hamiltonians that admit local errors $F$ with this property are very special. In particular, this assumption would not be satisfied if we tried to simulate single-site Pauli errors and take the energies as given by the Hamiltonian (4). However, we will work in the perturbative limit, where the effective Hamiltonian (29) is a good approximation to the true dynamics, up to a close-to-trivial unitary transformation. When restricted to the lowest energy eigenspace of the unperturbed Hamiltonian, this effective Hamiltonian reduces to (30), a stabilizer Hamiltonian, where the spin-boson simulation methods can easily be applied for local Pauli errors. As such, we will choose errors $F$ to be arbitrary single-qubit Pauli operators. The main subtlety is that these “single-site” errors we apply in the effective model actually correspond to quasi-local operators in the physical model, but this does not seem physically unreasonable, and in any case strict locality is restored in the perturbative limit.

By restricting to the lowest energy eigenspace of the unperturbed Hamiltonian, we would be significantly restricting the possible errors that may occur. In some of the regimes we consider, the system will (almost) never leave this eigenspace over the timescales of interest. However, it is also desirable to provide some account of the dynamics of our system outside these special regimes. In order to do this, we note that in a system where the energy spectrum is suitably close to that of a stabilizer Hamiltonian, we might nevertheless expect to approximate the thermalization process on small lattices over relatively short timescales and low energy scales by similar Monte-Carlo methods.

In particular, if an error $F$ takes Hamiltonian eigenstates to states whose spectral support is mostly contained within a small energy band, we can use the energy expectation value as an approximation to calculate $\Delta_{F}$. On general grounds, we expect that in the perturbative limit and at low energies, local Pauli operators $F$ will jump between narrow spectral bands in our system. A simulation where $\Delta_{F}$ is approximated by the difference of energy expectation values will then be reasonable so long as the total number of jumps is also small (i.e. small timescales and lattice sizes), though we leave a quantitative justification of the accuracy of such a simulation for future work.

IV.2 Error-correction protocols

In standard topological error correction protocols, the system is imagined to be initialized in a ground state of the Hamiltonian, before undergoing some noise channel, often a model for thermal noise. Local measurements are made on the system to determine an error syndrome, which is then fed into a decoding algorithm that will determine an appropriate recovery operation to return the system to the code space. Our protocol is identical to this setup, except for two small subtleties. The first is that the pseudo-ground space of our model is used as a code space in lieu of a truly degenerate ground space. The second is that we imagine measuring local operators in the effective model of Eq. (29). These operators correspond to quasi-local operators in the original Hamiltonian (4), as opposed to strictly local measurements.

Our decoding algorithm is based on the standard Perfect Matching Algorithm (PMA) decoders for the toric code Edmonds (1965); Kolmogorov (2009). These can be applied with only minor modification in our setting. Recall that in the limit $J_{x},J_{y}\ll J_{z}$, the KHM approaches the toric code Hamiltonian, where each toric code qubit corresponds to the 2-dimensional $K^{z}_{\mathbf{q}}=+1$ eigenspace at a given dimer of the honeycomb model. Before being able to apply a standard toric code decoding algorithm, our decoder must return the system to the common $K^{z}_{\mathbf{q}}=+1$ eigenspace at each edge.

Explicitly, the error correction protocol first takes measurements of all $K^{z}_{\mathbf{q}}$ and $W_{\mathbf{q}}$ operators. Following these measurements, any $K^{z}_{\mathbf{q}}$ operator that gave a $-1$ measurement outcome is returned to its $+1$ eigenspace by the application of an $X$ operator on either of the qubits of the corresponding edge (for concreteness, we choose the upper qubit). This will of course modify the measured $W$ syndrome so that it is no longer an accurate representation of the state, though the updated $W$ syndrome can easily be calculated by the Gottesman-Knill method Gottesman (1999) without needing to remeasure the system. Within the $+1$ eigenspace of all $K^{z}_{\mathbf{q}}$ operators, the $W_{\mathbf{q}}$ operators correspond to the stabilizers of the toric code, and so the updated $W_{\mathbf{q}}$ syndrome is passed to the PMA toric code decoder which returns an appropriate final recovery operation. Since the details of a PMA decoder for the toric code are well studied, we will not describe this algorithm, and refer the interested reader to e.g. Ref. Dennis et al. (2002) and references therein. Following this procedure, we are guaranteed to return to a state in the common $+1$ eigenspace of all $K^{z}_{\mathbf{q}}$ and $W_{\mathbf{q}}$ operators. Notice that this is precisely the ground space of the effective Hamiltonian (29). We regard the error correction as having been successful if the combined noise and recovery operations commute with the logical operators of the code. While this is far from an optimal decoding strategy, its effectiveness is sufficient for this exploratory study, and its close relation to the analogous toric code protocol makes it suitable for direct comparison of results.

One final remark is that in these simulations, we neglect the coherent sources of error studied in Sec. III, such as dephasing due to inexact degeneracy of the ground space. Since we have proved that these effects are exponentially suppressed in the system size, neglecting these additional errors is justified for our desired precision.

IV.3 Numerical results

We conduct simulations over a broad range of temperatures, focusing on several distinct regimes. The first regime we test is the ultra-high-temperature limit, in which the rates $\gamma_{F}\approx\frac{1}{\beta}$ are independent of the error $F$ and so we may effectively disregard the Hamiltonian. This serves as a benchmark of our model, corresponding to the well studied case of i.i.d. depolarizing noise on each spin.

Following this benchmark, we consider three temperature regimes: low, intermediate and high. The low temperature regime for our system corresponds roughly to $\beta\gtrsim 10^{5}$, where we find interesting finite-size effects beginning to become significant. The intermediate temperature regime corresponds roughly to $\beta\sim 10^{4}$, and is chosen such that $\beta\Delta_{F}\sim 1$ for errors $F$ that create two vortices. In the low- and intermediate-temperature regimes, for the lattice sizes studied and over timescales comparable to the memory lifetime of the system, the probability of encountering a broken dimer is negligible. The high temperature regime has $\beta\sim 10^{1}$, which is chosen such that $\beta\Delta_{F}\sim 1$ for errors $F$ that create a broken dimer. At these temperatures, the density of $-1$ $K^{z}_{\mathbf{q}}$ eigenvalues becomes measurable during our simulations, and so the mean-energy approximation is required.

All numerical results presented in this work have parameters $J_{z}=\frac{1}{2}$ and $J_{x}=J_{y}=0.1$, and energies are taken by the effective Hamiltonian (29) calculated to $6^{\rm th}$ order. For convenience, we have also taken different boundary conditions to those used in the previous sections. Instead of the lattice of Fig. 1, we will now use the rotated boundary conditions shown in Fig. 2, and denote the corresponding linear lattice size by $\tilde{N}$.

IV.3.1 The ultra-high temperature benchmark

In the infinite temperature limit, the Hamiltonian plays no role in the thermalization dynamics, and so we can approximate these dynamics by an i.i.d. depolarizing channel for each qubit. We can also neglect coherent evolution in this regime, since it is assumed to act on a much longer timescale than the memory lifetime, which will be extremely short. As such, this can be treated as a benchmark for the honeycomb model and readily compared to error correction of the standard toric code under i.i.d. depolarizing noise.

The results of this simulation are shown in Fig. 3, where the threshold value is given as a memory lifetime instead of the typical depolarizing channel strength, for more accurate comparison with later results. A critical memory lifetime (defined as the time below which increasing the lattice size reduces the error rate) can be found at around $\frac{t_{c}}{\beta}\approx 0.04$. Noise and error correction for the honeycomb model in the ultra-high temperature regime are extremely similar to the standard toric code, since the existence of the honeycomb Hamiltonian plays no role except to define the initial encoded states. The only subtlety is that the depolarizing channel applied to the honeycomb model does not correspond to the depolarizing channel applied to the toric code, as toric code qubit operators are encoded in two qubits of the honeycomb model. As such, since we sample single-site Pauli operators independently, we effectively consider a noise source for the toric code that biases towards phase-flip errors over bit-flip errors.

On small timescales (equivalently small error rates), the simulation time is approximately related to the depolarizing channel strength as $p_{\mathrm{dep}}\approx\frac{3}{\beta}t$ (since there are three possible Pauli error events per qubit, each with rate $\frac{1}{\beta}$). This allows us to compare our memory lifetime to the standard toric code under depolarizing noise using a PMA decoder, which has a threshold of around $p_{\mathrm{dep}}^{*}\approx 0.158$ Dennis et al. (2002), corresponding to a memory time of $\frac{t_{c}}{\beta}\approx 0.05$. Note also in this comparison that since each of the physical qubits of the toric code corresponds to two physical qubits in the honeycomb model, the effective error rate per encoded toric code qubit is approximately double that of the honeycomb qubits.

IV.3.2 Low temperatures

For low temperatures, $\beta\gtrsim 10^{5}$, we find that the appearance of broken dimers over the memory lifetime of the system (for the lattice sizes we consider) is so infrequent that they were never observed. The error rate as a function of time is shown for various lattice sizes in Fig. 4.

We do not find evidence of threshold behavior in this regime. In fact, it appears that the pseudo-threshold, representing the time below which increasing the size of the system reduces the logical error rate, is itself increasing with system size. This can be understood as a finite-size effect in this low temperature regime. It is known that at low temperatures and on finite size lattices, there is a regime where the random walks of a single pair of toric code excitations (vortices in the honeycomb model) are the dominant source of error Freeman et al. (2014); Brown et al. (2016). This effect leads to increasing suppression of the logical error rate as the lattice size decreases, as seen in Fig. 4.

At such low temperatures, since it is unlikely that any error would lead to a broken dimer, and since we assume a local coupling of each spin to the environment, this effectively restricts the errors considered in the honeycomb model to Pauli $Z$ errors (all others occur with vanishing probability). It is easy to see that these errors can only move vortices horizontally along 1-dimensional rows, as opposed to the typical situation for the toric code where the excitations may be moved in any direction. As such, we can intuit the dynamics of excitations by considering them to be undergoing a 1-dimensional random walk, instead of the 2-dimensional walks seen in Refs. Freeman et al. (2014); Brown et al. (2016).

As in the conventional toric code, this low temperature regime is expected to persist while the number of anyons is low, i.e. $N^{2}\mathrm{e}^{-\beta\delta E}\lesssim 1$ for $\delta E$ the gap, giving a critical $\beta^{*}\sim\frac{\ln N}{\delta E}$. Within this regime, the memory lifetime has been argued phenomenologically to be inversely proportional to the probability of a topologically non-trivial walk, $\Pi_{\mathrm{nt}}$ Freeman et al. (2014). In the 2-dimensional case, numerical evidence Brown et al. (2016) supports this scaling with $\Pi_{\mathrm{nt}}^{2D}\sim\frac{1}{\ln N}$, while for our 1-dimensional walk behavior, we expect $\Pi_{\mathrm{nt}}^{1D}\sim\frac{1}{N}$ Freeman et al. (2014). This would increase the memory lifetime by a factor of $\frac{N}{\ln N}$ in this regime compared to the behavior of the standard toric code, strengthening the effect of this error-suppression mechanism on larger lattices. This interesting finite-size effect in the perturbative regime of the honeycomb model may warrant further study in order to determine under what circumstances it could be exploited to produce improved memory lifetimes or logical error rates. If this were desirable, the decoding algorithm could presumably also be optimized to take this extra structure into account.

This effective 1-dimensional walk behavior is an artefact of our having restricted to independent errors on each qubit. Typically this is justified by the expectation that physical errors will act quasi-locally, and that there will be little qualitative distinction between strictly local and quasi-local error models. However, it is clear that in the present case, if we were to allow suitable 2-local errors, we could engineer behavior equivalent to the standard toric code under depolarizing noise. Thus a more careful consideration of the physically relevant noise model for a given implementation must be made when considering error correction in the low-temperature regime of the honeycomb model.

IV.3.3 Intermediate temperatures

At intermediate temperatures ($\beta\sim 10^{4}$), we are still in a regime where broken dimers are sufficiently unlikely that they are unobservable for the system sizes and timescales we simulate. However, it differs from the low-temperature regime in that the failure modes are no longer dominated by randomly walking anyons. Rather, they are dominated by the more typical situation where clusters of errors percolate beyond the point at which they can be reliably distinguished.

The results of our simulations are shown in Fig. 5. The memory lifetimes (or equivalently error thresholds) are estimated as the common crossing point of the curves for different lattice sizes, and plotted in Fig. 6. At lower temperature and small lattice sizes (most notably Fig. 5d at $N=16$), we see deviation from this common crossing point, as the finite-size effects seen in the low-temperature regime (Fig. 4) begin to become significant. Nonetheless, at larger lattice sizes the expected threshold behavior is observed. The memory lifetimes are found to be exponential in inverse temperature, as is the case for the standard toric code Brown et al. (2016).

IV.3.4 High temperatures

The low- and intermediate-temperature regimes considered in the previous sections were characterized by the fact that broken dimers were unlikely to appear on timescales shorter than the memory lifetime. This is clearly unrealistic as we approach the thermodynamic limit, so in order to characterize the impact of leaving this regime, we consider two sets of effects—coherent diffusion of broken dimers, and a heuristic simulation of thermalization processes.

Dimer diffusion

We consider a state with a single broken dimer ($-1$ eigenstate of $K^{z}_{\mathbf{q}}$) and no vortices. We then simulate the time-evolution of this state under the effective Hamiltonian on dimer degrees of freedom (taken to $2^{\rm nd}$ order), before calculating the distribution of broken dimer positions, and hence the average distance travelled. This will allow us to determine whether the diffusion of broken dimers is significant on timescales relevant for error-correction. The results are shown in Fig. 7.

Mean-energy analysis

Finally, we perform a heuristic simulation of thermalizing dynamics at high temperatures, in which we approximate the quantum jumps of the monte-carlo method described in Sec. IV.1 by single-qubit Pauli errors with transition probabilities determined by mean energy differences in the effective Hamiltonian. This approach is expected to be accurate in the perturbative regime for small lattices at short timescales. The results are shown in Fig. 8.

As in the intermediate temperature regime, we see clear evidence of threshold behavior in the memory lifetime, though the computational cost of these simulations is higher and so we are restricted to probing slightly smaller lattices. The memory lifetimes are plotted in Fig. 9, where we see behavior consistent with linear lifetime in inverse temperature, in contrast to the intermediate temperature regime. This is a result of the small $\beta$ limit of Eq. (140) giving the error rate for vortex-creation errors as $\gamma\sim\frac{1}{\beta}$. Note however that the dimer-breaking error rates are not yet in this high-temperature linear-scaling regime. This leads to a contrast between the lifetimes in the ultra-high temperature regime (Fig. 3) where all errors are in the linear-scaling regime, and $\frac{t_{c}}{\beta}\approx 0.04$, and the temperatures considered in Fig. 8, where only Pauli $Z$ errors are in the linear-scaling regime, and correspondingly the lifetimes are increased to $\frac{t_{c}}{\beta}\approx 0.14$.

We also plot the observed density of broken dimers over relevant timescales in Fig. 10. We find dimer density equilibrating exponentially fast at high temperatures, as expected. The density of dimers is always below $10\%$ even in our highest temperature simulations, and for $\beta\geq 20$, we observe no broken dimers for lattice sizes and timescales studied.

It should also be noted that the diffusion behavior studied in Fig. 7 is unlikely to significantly affect the error correction properties of this model, since at the higher temperatures where the dimer density is non-negligible (e.g. Fig. 10a), the memory lifetime is so short that the diffusion distance is negligible, and at lower temperatures where the memory lifetime is long enough to give average diffusion distance on the order of $0.1$ lattice spacings (e.g. Fig. 10b), the density of broken dimers is so low that their diffusion is unlikely to cause them to interact or play a major role in logical failures.

To summarize, we do not find evidence that the appearance of broken dimers at high temperature has any qualitative effect on the honeycomb code on relevant timescales, and we observe behavior similar to that of the standard toric code.

V Discussion

Our analysis and numerical simulations highlight behaviors in the honeycomb model that differ in several ways from error correction in the standard toric code. Despite this, our results suggest that these disparities will not qualitatively affect the properties of the honeycomb model as a quantum memory.

Of course, as shown broadly in Sec. III.1, the key error correction properties of degeneracy and local indistinguishability hold only approximately for the honeycomb model. Analytically, we have provided an explicit exponential bound on the energy splitting of the codespace, which is important for bounding the dephasing time of the memory, as well as on its local indistinguishability and approximate correctability, which are crucial for determining the parameters of the code. By taking advantage of the special structure of the honeycomb model, these bounds are sharper that what appears to be possible using the known toolkit for general topological phases, and have the added bonus of explicit constants.

Numerically, we have found an exponential lifetime in inverse temperature, as holds for the toric code. We also found evidence that the appearance of non-topological broken dimer excitations has little qualitative effect on the error correction process. In regimes where the memory lifetime would be long enough for dimer diffusion to be significant, the temperature is so low that dimers are exceedingly unlikely to arise. While this is quite fortuitous, it is unclear to what extent this particular phenomenon would extend out of the perturbative regime, where our simulation methods are unable to probe. In any case, we expect that error correction in general will be possible everywhere in the gapped phase, though of course the memory lifetime is expected to decrease as the gap decreases. How to reasonably implement the detection and manipulation of anyons has been considered in various experimental contexts, most thoroughly for optical lattices Zhang et al. (2007); Pachos (2007); Jiang et al. (2008); Aguado et al. (2008); Dusuel et al. (2008).

As in the standard toric code, we identified several distinct regimes of behavior, the most interesting of which is the low-temperature regime, where finite size effects dominate the error modes. While this kind of regime exists in the standard toric code, we argued that the structure of the honeycomb model leads to an improved lifetime scaling by a polynomial factor under local thermal noise. Again, while fortuitous, this effect relies on the perturbative structure of the honeycomb mode, and so it is unclear how far outside the perturbative regime it can be expected to persist. This may warrant further study to determine whether this effect can be exploited in near-term experimental implementations.

While these departures from the behavior of the toric code are important to understand when designing good error correction protocols, none of these novel behaviors are particularly shocking, nor would they seem to be particular obstacles to successful error correction in the honeycomb model.

Studies of approximate topological quantum error correcting codes and their behavior are in their infancy, as are designing tools to probe these systems, and developing protocols to ideally compensate for the deviations from exact codes. There are several issues that must be considered when analyzing such systems, one of which is the motivating experimental setup.

In typical studies of topological error correction in idealized models, the error correction protocols are designed to return precisely to the ground space of the Hamiltonian, as we have considered above. A clear motivation for such a protocol is that the system could be initialized by cooling to a state exactly within the ground space, and the error correction procedure is fine-tuned to return to a state within this space. In a more realistic system, any imprecision in the implementation or knowledge of the Hamiltonian, or any restriction on the ability to perform non-local unitary operations perfectly will prevent a return to the ground space (i.e. the overlap between the recovery space and the Hamiltonian ground space will generically decay exponentially quickly in the system size). One might imagine that after performing a recovery operation that returns to a state “close” to the ground space, the experimentalist can again cool to an exact ground space state, but for 2D topologically ordered systems it is not clear that this can be done without performing a logical error.

As such, when considering topological error correction in realistic systems, there are several kinds of error correction protocols that may be of interest. The first is as described above, where the system is initialized into the ground space, and the error correction algorithm returns to this space. The second involves again initializing the memory into the ground space, whereas the error correction protocol attempts to return not to this initial space, but to another code space. A reasonable candidate for this space in some circumstances is the ground space of an RG fixed point of the same phase. The third class of error correction protocols is more suited to multi-round error correction, and does not initialize in the ground space of the Hamiltonian. Instead it initializes in the same space in which it intends to return (such as the RG fixed-point ground space).

It may be of interest to consider the properties of protocols falling within these other error correction strategies, depending on the physical application in mind. Furthermore, if the intended setup has a code space other than the ground space, then it is not only of interest to understand the degeneracy splitting and local indistinguishability of the ground space, but also the dephasing rates and code properties of this other space. A fundamental question to consider is then the role of the Hamiltonian in these kinds of topological quantum information storage protocols. Although the Hamiltonian would normally provide some level of passive error suppression through the presence of a mass gap, this is no longer obviously the case if we do not always encode information in the ground space of the model. Nonetheless, this is presumably the practically relevant scenario for an experimentalist attempting to use a real topologically ordered system for information storage in the laboratory.

An obvious extension of this work is to consider the honeycomb model as a quantum memory in the full fault-tolerant setting, where faulty measurements are allowed. It would also be very interesting to simulate the model outside the perturbative regime and/or in the non-abelian phase of the model. Though it is not obvious how to perform such simulations efficiently, recent progress in tensor network descriptions of the honeycomb model Schmoll and Orús (2017) and tensor network simulations of noise and error correction Darmawan and Poulin (2017) may be a fruitful avenue for such research.

We believe that many of the results in this paper may give leading clues to behavior in more general approximate topological codes. The structure of the honeycomb model has provided us a starting point to investigate more realistic quantum memories than the exact toric code, and developing methods to analyze more general systems is an important step in understanding topological quantum error correction in realistic devices.

Acknowledgements.

We thank Jens Eisert, Min-Hsiu Hsieh, Tobias Osborne, David Poulin, and Ted Yoder for valuable discussions. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. CGB was supported by the ERC grant QFTCMPS and by the cluster of excellence EXC 201 Quantum Engineering and Space-Time Research. STF was supported by the Australian Research Council via EQuS project number CE11001013 and by an Australian Research Council Future Fellowship FT130101744. This is an author-created, un-copyedited version of an article published in the Journal of Statistical Mechanics: Theory and Experiment. IOP Publishing Ltd is not responsible for any errors or omissions in this version of the manuscript or any version derived from it. The Version of Record is available online at doi.org/10.1088/1742-5468/aa7ee2.