Exact description of the boundary theory of the Kitaev Toric Code with open boundary conditions
Yevheniia Cheipesh Lorenzo Cevolani Stefan Kehrein






Institut für Theoretische Physik,
 Georg-August-Universität Göttingen, 37077 Göttingen, Germany




 Institut für Theoretische Physik,
 Georg-August-Universität Göttingen, 37077 Göttingen, Germany






Institut für Theoretische Physik,
 Georg-August-Universität Göttingen, 37077 Göttingen, Germany






In this work we consider the Kitaev Toric Code with specific open boundary conditions. Such a physical system has a highly degenerate ground state determined by the degrees of freedom localised at the boundaries. We can write down an explicit expression for the ground state of this model. Based on this, the entanglement properties of the model are studied for two types of bipartition: one, where the subsystem A is completely contained in B; and the second, where the boundary of the system is shared between A and B. In the former configuration, the entanglement entropy is the same as for the periodic boundary condition case, which means that the bulk is completely decoupled from the boundary on distances larger than the correlation length. In the latter, deviations from the torus configuration appear due to the edge states and lead to an increase of the entropy. 
We then determine an effective theory for the boundary of the system. In the case where we apply a small magnetic field as a perturbation the degrees of freedom on the boundary acquire a dispersion relation. The system can there be described by a Hamiltonian of the Ising type with a generic spin-exchange term.


Exact description of the boundary theory of the Kitaev Toric Code with open boundary conditions
    Stefan Kehrein
    December 30, 2023
===============================================================================================





§ INTRODUCTION


The Landau theory of phase transitions is one of the most precise paradigms in modern physics. Once a system undergoes a symmetry breaking, the macroscopic properties of its ground state change.  Some phase transitions are not characterised by the local order parameter, instead the topological order can be used <cit.>. The latter plays a crucial role in the existence of so–called topological insulators. In these systems, conduction can take place only along the boundary, while in the bulk one has an insulating state. Topological properties are extremely robust, they cannot be destroyed by simple perturbations. 
Many of these states, e.g. the Fractional Hall Effect, involve the correct description of degrees of freedom that are not located in the bulk but on the boundaries of our system. If we consider for example a bi–dimensional system with open boundary conditions, a large number of degenerate ground states will be located on the boundary. These states are called edge states <cit.> and their correct description can be complicated. Physically, these states are responsible for the behavior of the system in the low energy regime. The bulk theory is usually gapped, while the boundary theory is known to be gapless, meaning that in the low energy limit the boundary states will affect the behavior of the observables. A remarkable example of such effect is the current quantization in Hall systems <cit.>.

In order to grasp the physics of such states, exactly solvable models as <cit.> play an important role. These can be used as playground to develop exact mathematical concepts that can then be extended to other more complicated systems. One of these models is the so–called Kitaev Toric Code (KTC) <cit.>, whose peculiarity is due to the toric geometry. This model can be used to store quantum information in its 4-fold degenerate ground states with different topological properties <cit.>.



In this paper we study the KTC with particular open boundary conditions: the coupling constants of some operators are set to zero along the boundary. The model is exactly solvable in this case.  
The exact representation of the ground state gives a clear expression of the edge states, which are localized on the boundary of the system, but also traverse the bulk. The specific form of the ground states can be constructed from different operators acting separately on these regions. The resulting state exhibits a clear connection between the bulk and the boundary. We then use the derived exact expression for the ground state to compute the bipartite entanglement entropy of our model. We find that this has the same value as for the periodic boundary condition case if the partition is in the bulk of the system. If the boundary is also partitioned, then the value of the entanglement entropy is increased by a factor that depends on the length of the boundary of a subsystem. Due to the fact that the geometry is topologically trivial, the ground states cannot be used to store information safely as in the original model. In the end, we study how the boundary degrees of freedom can be excited using a magnetic field to induce a non-flat dispersion relation.

The paper is structured as follows: in the Section I we describe the model and the boundary conditions we choose, in the Section II we focus on the construction of the ground state and on its interpretation on the bulk–to–boundary correspondence, the third Section is dedicated to the computation of the entanglement entropy, and in the last Section we focus on the excitations of the boundary theory and their dispersion relation.




§ THE KITAEV TORIC CODE WITH OPEN BOUNDARY CONDITIONS


We want to study a two–dimensional spin system with open boundary conditions (OBC). The system is a m× n rectangle with spin-1/2 located at the links between the lattice sites which Hamiltonian reads as:

    Ĥ   = -J_e∑_vÂ_s - J_m∑_pB̂_p, 
    Â_s    = ∏_starσ̂^x_i,     B̂_p = ∏_plaquetteσ̂^z_i.

B̂_p is called "plaquette" operator, and Â_s is the "star" operator. These two operators trivially commute and can have eigenvalues ± 1. The ground state of this model is an eigenstate for all the star and plaquette operators with eigenvalues +1.  The Hamiltonian (<ref>) with periodic boundary conditions (PBC) is called the Kitaev Toric Code (KTC) and has already been studied in several works in the literature <cit.>. Its ground state is 4-fold degenerate and has non–trivial topological properties living on a torus due to the PBC. This model has been proposed to be used as an error–free memory <cit.>. Various more complicated models with the same properties have since been discussed in the literature <cit.>.


The specific case of the Hamiltonian (<ref>) with different OBC has been studied in several cases <cit.>. These models are referred to
as planar codes. One of the main results of these works is the determination of the edge states which are all degenerate.

In our work we are going to consider a modified KTC with particular OBC: we set to zero the J_e constants along the boundary of our system, see Fig. <ref> (B) for a graphical representation. This means that the plaquette operators are still active on the boundary, while the star operators are not. We will discuss the physical implications of turning on and off these operators in the next section. Clearly, it is also possible to set the plaquette operators to zero instead, see Fig. <ref> (A). From now on we will always refer to the configuration (B), but our results can be extended to the other case by using the eigenvectors of σ_x as a local basis instead of ones of σ_z as we will do henceforth.




 
 





 



§ GROUND STATE PROPERTIES OF THE KTC WITH OPEN BOUNDARY CONDITIONS

 Let us first take a look at the ground state of the original KTC <cit.>. The full set of degenerate ground states can be constructed starting from the trivial state when all the spins are down


    |0⟩=∏_j|↓⟩_j

where j spans the full lattice.

Then the ground state of the KTC in presence of PBC has been demonstrated to be <cit.>:

    |GS⟩_a = √(2)∏_j( 𝕀 + Â_j/√(2))Q̂_a|0⟩

where the index j runs over the different lattice sites. The operators Â_j were defined in Eq. (<ref>). The operator Q̂_a belongs to a set of 4 different operators {𝕀, Ŝ_1, Ŝ_2, Ŝ_1 Ŝ_2 }, where Ŝ_1, Ŝ_2 correspond to closed loops that flip spins along two main circles of the torus, correspondingly. One sees that the ground state with PBC is four–fold degenerate.

In the case of OBC, the previous state is still a valid ansatz for the ground state. However, it is just one among many degenerate states that is also valid in the case with PBC. The degrees of freedom obtained by removing the constraints on the boundary enrich the ground state subspace  compared to the one  in the original model. 


Edge states arise on the boundary creating a huge degeneracy and can be written as:

    |GS⟩ = ∏_j^Σ(𝕀+ Â_j/√(2))∑_e⃗ c_e⃗∏_i^ℒ-1( 𝒲_i )^e_i|0⟩.

The operators 𝒲_i= σ_i^xσ_i+1^xσ_α^x are represented in Fig. <ref> and they involve just the degrees of freedom located at the boundary. 

The vector e⃗ has length equal to the boundary length ℒ-1 and every entry can be either 0, which means no operator 𝒲_i on that point,  or 1, which means that 𝒲_i is "activated" on that point. One should note, however, that here we consider rectangular lattice which means that we have additionally 4 corners that have a different configuration compared to the sites along the boundary. That is why, strictly speaking, we should also add 4 additional operators 𝒲_j = σ_j^xσ_j+1^x that  have the same nature and structure as ones introduces for the boundary and represented in Fig. <ref>, but involve only 2 spins that are  flipped on the corner. This is necessary in order to full fill the condition that the product of all  𝒲_i is trivial. In the case of circular geometry one does not have this issue. The number of ground states is then equivalent to the number of different vectors present: 2^ℒ-1. The coefficients c_e⃗ take into account all the possible superpositions of these boundary states. It is important to notice how the state, previously written in Eq. (<ref>), takes the form of an expansion of states that are factorised into a part defined on the boundary, the one involving 𝒲_i operators, and a part defined on the bulk, the one involving Â_i operators.

Eq. (<ref>) is one of our main results, in the following two sections we will describe how to obtain this state starting from the operators defining Hamiltonian (<ref>) in full details. The reader not interested in these details can skip this part and go directly to the Section IV dedicated to the entanglement entropy.




 §.§ Construction of the ground state


Before we focus on the structure of the ground state, we want to demonstrate that it is massively degenerate. The degeneracy of the ground state can be determined by counting the degrees of freedom of the system and then subtracting the number of constraints imposed by the stabilizers <cit.>.
For the sake of simplicity, let us consider an  N× N lattice. There, the number of q-bits is 2N^2 + 2N, where  2N comes from the boundary. The number of plaquette operator is N^2 while the one of the stars is N^2 - 2N + 1, due to the fact that they are absent on the boundary. So the total number of all configurations allowed by the constraints of the ground state is 2^4N-1=2^ℒ-1 where ℒ=4N is the length of the boundary.

Now we want to construct the ground state itself following the requirement that it should be an eigenvector of all Â_s, B̂_p. Such a vector exists as far as [ Â_s,B̂_p] = 0. We use the local basis of σ_z, where

    |ϕ_1⟩ = ( [ 0; 1 ])     |ϕ_2⟩ = ([ 1; 0 ])



    σ̂_x|ϕ_1⟩   =|ϕ_-1⟩    σ̂_x|ϕ_-1⟩=|ϕ_1⟩
    σ̂_z|ϕ_1⟩   =|ϕ_1⟩    σ̂_z|ϕ_-1⟩=-|ϕ_-1⟩.

This is the most suitable choice of local basis for our boundary conditions.

Let us now proceed to build the ground state. To do so, we need to have a clear picture of what the star and plaquette operators physically represent.

It is useful to start from considering B_p, since it is diagonal in our local basis. For the unit cell with four spins as in Fig. <ref>, the product of the spins along the square has to be equal to +1. Clearly, only the configurations with an even number of spins up satisfy the ground state condition. This requirement  will be important also for the construction on the boundary since the plaquette operators are the only ones active there. More interesting is the action of the star operators. Since each star operator is composed of four σ_x matrices, it corresponds to a simultaneous flip of all the corresponding spins. One can see the action of star operators on all possible states in Fig. <ref>. In other words, states are transformed into each other under the action of the start operator. It is important to note, that the new state is still an eigenstate of a plaquette operator with an eigenvalue +1 since the star operator cannot change the parity of spins.



The previous considerations can be used to construct the requirements for the ground state of the KTC (namely that each plaquette can have only an even amount of spins flipped) and also to understand the action of the star operators on such states.

Now we explicitly construct the ground state of the model with OBC. Let us start considering the boundary of our system where only plaquette operators act. Starting from a configuration where all the spins are in the eigenstate +1, we can see that the path of spin flips connecting two points on the boundary does not cost any energy as we flip two spins for every plaquette, see Fig. <ref> (A).  It means that the ground state consists not only of the state |0⟩ where all the spins look up but can also have "paths" of flipped spins that end on the boundary.



Since the star operators transform states into each other as in Fig. <ref>, the ground state should include all the possible states, allowed by plaquette conditions, with equal coefficients. In this case we stay in the same orbit of the group, generated by all the possible products of the star operators. As far as this group acts only in the bulk of the system, because of our particular boundary conditions, it does not change the spin configuration on the boundary which means that the element of the basis |l⟩ from which the ground state is composed is defined by the particular set of points on the boundary. Thus we can build the ground state as:

    |GS⟩=∑_l c_l ∏_i ( 1̂ + Â_i ) |l⟩.

An important physical interpretation of the star operators can be seen by comparing Fig. <ref> (A) and (B): applying the star operator to the part of the path that connects the two boundaries changes the state, but just in the bulk. This connection between the Â_i operators and the modification of the paths (stretching) can be used to give a different form of Eq. (<ref>).

The structure of the deformations of the paths induced by Â_i can be recast in a more formal way <cit.>. We introduce the group 𝒜 that contains all the possible stretches of any given path |l⟩. The elements of 𝒜 are operators ĝ generated by all possible products of the star operators in the bulk of the system. The ground state can then be rewritten as a superposition of the action of all the elements of this group:

    |GS⟩ = 1/√(|𝒜|)∑_l a_l ∑_𝒜ĝ|l⟩,

where we used the fact that any configuration of the star operators is contained in 𝒜, meaning that ĝ∈𝒜→Â_i ĝ = ĝ^', and we normalised the vector using |𝒜|=2^Σ where Σ is the total number of independent star operators. Eq. (<ref>) is even easier to interpret physically than Eq. (<ref>): the boundary and the bulk are completely decoupled due to the absence of stars that can stretch spin flips from the bulk to the boundary and vice versa.


Eq. (<ref>) says that if we have two points on the boundary, we can connect them and, in order to be a ground state, the resulting state should be a superposition of all the possible "stretches" of the paths between these two points. We can identify them as one equivalence class, namely, we say that two configurations are equivalent if they correspond to the same set of points on the boundary. This allows us to work only with one element of each class and detaches ourselves from the consideration of the bulk.

This particular set of paths is represented in Fig. <ref>, where we simply connect the two points on the bulk by the path closest to the boundary. The action of the elements of 𝒜 will take care of modifications in the bulk to have a proper ground state. 

We want to focus now on the boundary states |l⟩. We will demonstrate that it is possible to create all of them starting from the simple building blocks. One of these building blocks is shown in Fig. <ref> which is represented by three spin flips connecting two nearest–neighbour spins on the boundary:













    𝒲_i = σ_i^x σ_i+1^x σ_α^x,

where i denotes the spins on the boundary and α the only one in the bulk.

The boundary of a two–dimensional system can be represented exactly by a one–dimensional system with periodic boundary conditions. All the physical quantities will depend just on the distance between i and j, which label the points on the boundary, and never on i and j alone.
These operators can be used to construct any type of boundary theory simply by multiplying them appropriately	. Let us take the example of Fig. <ref> where a segment connecting two points is created:

    𝒲_i( R=3 ) = 𝒲_i𝒲_i+1𝒲_i+2.





We can check that all the ground states are built by these operators. If the boundary length is ℒ, we have ℒ 𝒲_i operators. Each of them can either be activated or not, creating different types of boundary states. The total number of possibilities is then 2^ℒ. Also, if we look at the boundary alone, the state where all the 𝒲_i are applied is the same to the identity, because:

    ∏^ℒ_i 𝒲_i= ∏^ℒ_i σ_i^xσ_i+1^xσ_α^x= 1̂_boundary(∏_Bulkσ_α^x).

This reduces the number of possible boundary states from 2^ℒ to 2^ℒ-1, which coincides with the number of ground states we obtained computing the constraints and degrees of freedom.


All the boundary states can then be written as functions of numbers e_i that are 0 if the corresponding 𝒲_i is not applied to |0⟩ and 1 if it is applied. This numbers can then be stored into a vector e⃗ which then represents the boundary theory as:

    |l⟩→e⃗  |l⟩= ∏_i^ℒ( 𝒲_i )^e_i|0⟩.

As we said before, the boundary conditions can restrict the dimension of these spaces. The boundary, when considered as a system by itself, has a pure periodic boundary conditions, meaning that the vector where all the 𝒲_i are present is the same to the one where none of the spins is flipped.

    (0,0,0, …,0,0) ∼( 1,1,1…,1,1 ),

which means that the effective number of degrees of freedom is ℒ-1 and the total dimension of the space is then 2^ℒ-1 which is exactly the total number of ground states of the model.


We can now focus on the coefficients a_l of the expansion in Eq. (<ref>). Since the boundary of our system is topologically equivalent to a one–dimensional system with PBC, a lot of these coefficients will be equal due to the translational symmetry. We can start from the fundamental building blocks, Eq. (<ref>) and Fig. <ref>. Since all these operators are equal, their coefficients have to be equal too: a_1=a.
The fact that all the other states in the superposition can be obtained by the multiplication of these fundamental blocks fixes the other coefficients. Since the multiplication of operators 𝒲_i  gives more complicated boundary state, the generic element of the expansion is given by

    |l⟩= ∏_i^ℒ( 𝒲_i )^e_i|0⟩  →   a_l = a^l e^iφ_l,

where l=∑_i e_i^2=|e⃗|^2, φ_l is the phase of the complex number, we fixed φ_1=0. Physically, this means that the states with the same number of operator 𝒲_i differ just by a phase.


We want to emphasise here, that even though the considerations made in this section were formulated for a square geometry with no star-operator on its boundary. Our conclusions can be easily extended to more exotic geometries and boundary conditions where some of the operators of the Hamiltonian are set to zero along some parts of the boundary. The ground state can be constructed starting from the 𝒲_i operators applied over a part of the boundary and by the application of the bulk operators. The dimension of the ground state depends on the length of the boundary where one has removed the stabiliser operators. 





§ ENTANGLEMENT ENTROPY OF THE GROUND STATE



Once we have the exact analytical expression for the ground state, we can now compute  the entanglement entropy of a bipartition of the system. The system is in a pure state, of which the density matrix is

    ρ̂=|GS⟩⟨GS|.

At zero temperature, the entanglement in a many–body system can be estimated using the von Neumann entropy <cit.>:

    S_A = -_A ρ̂_A lnρ̂_A,

where ρ̂_A=_Bρ̂. The entanglement entropy of the KTC has been computed exactly in <cit.>

    S_A ≡S̃_A = (ℒ_A/B-1)ln 2.

These calculations show that the entanglement entropy follows the area law, as expected from the presence of a finite correlation length <cit.>, plus a constant term which can be computed exactly. This constant term is the topological entanglement entropy <cit.>, and it is connected to the quantum dimensions of the anyons present in the theory.

We want to compute the ground state entanglement entropy for our system with the specific boundary conditions specified in Fig. <ref>. The rectangular lattice is split into two sub–systems (A and B) as presented on the left panel of Fig. <ref> so that every spin lies either in sub–system A or B, no spins are shared between them. With respect to the KTC, we have more freedom in the choice of the partition between the subsystems A and B. In the particular case when a subsystem A is smaller than the subsystem B, we can consider two main configurations: the first is where the boundary of the system is completely contained in the subsystem B, and the second is when the two subsystems share part of the boundary. 

The first case is similar to the case of the KTC with PBC. There, the information about the boundary states is completely destroyed by the trace over the subsystem B (we are assuming that A is at least one correlation length away from the boundary). This and the fact that the bulk and the boundary are completely decoupled due to our boundary conditions makes it intuitive that the entanglement entropy takes the same form as for PBC, S_A=S̃_A.

However, when  also the boundary is partitioned, we have an increase in the entanglement entropy due to the edge states:

    S_A = S̃_A + ln1 f(A)f(B),

where the function f(x) is a positive monotonously increasing function that depends only on the length of the boundary in the subsystem.

These are the main results of the section.The reader interested in the detailed derivation is encouraged to follow the next subsection. In other case, proceed to section V.


 §.§ Entanglement entropy of the boundary

We now want to focus on the second case, where two subsystems share part of the boundary, see Fig. <ref>. In this case, the boundary states will not be completely traced out performing _B, but they will leave their fingerprint on the entanglement entropy.


We start representing the ground state using the group 𝒜:


    |GS⟩ = 1/√(|𝒜|)∑_la^n_le^iφ_l∑_g ∈𝒜g  ∏_i^ℒ( 𝒲_i )^e_i|0⟩.


Let's first calculate the reduced density matrix. Before doing this, we should specify carefully what happens to the stars that are shared by two subsystems A and B (see Fig. <ref>). We decided to place the boundary of the two subsystems in between the spin and the lattice site. In doing this we cut some of the star operators between the two subsystems, see Fig. <ref>. Once we trace out the degrees of freedom of B, this operation fixes part of the spins in A, due to the ground state conditions. Keeping track of this reduction is fundamental to take properly into account the degrees of freedom in between A and B, which is then responsible of the entanglement entropy.

We start by tracing out the degrees of freedom of the subsystem B to obtain the reduced density matrix of the system A:


    ρ̂_A = 1/|𝒜|∑_l,ka^n_la^n_ke^i(φ_l - φ_k)∑_g,g' ∈𝒜_B (⟨0|_B ∏_i^ℒ( 𝒲_i )^e_i g_Bg_A/B^sharedg_A|0⟩_A ⟨0|_Ag'_Ag'_A/B^sharedg'_B  ∏_i^ℒ( 𝒲_i )^f_i|0⟩_B).


The previous expression can be simplified using two important facts: The boundary states in B subsystem are orthogonal, and the bulk states are orthogonal if their modifications are different, meaning:


    ⟨0|∏_i^ℒ( 𝒲_i )^e_i∏_i^ℒ( 𝒲_i )^f_i|0⟩ = ∏_i δ_e_i,f_i
       ⟨0| g_B^' g_B |0⟩ =δ( g_B^' -g_B ) .

Using these identities the previous equation can be simplified to

    ρ̂_A =
    1/|𝒜|∑_l,k∑_g_B,g_B^'∈𝒜_Bδ^n^B_l_n^B_kδ(g_B-g'_B) a^n_la^n_ke^i(φ_l - φ_k)
       ∑_g_A, g_A^'∈𝒜_Ag_A∏^A/Bσ̂^x ∏_i^ℒ_𝒜( 𝒲_i )^e_i|0⟩_A⟨0|_A ∏_i^ℒ_𝒜( 𝒲_i )^f_i∏^ A/Bσ̂^xg^'_A.

Where we used the fact that |0⟩=|0⟩_A⊗|0⟩_B because it is a product state.

We can then divide the full group of modification of the paths in two different subsets 𝒜_A and 𝒜_B containing all the stars fully contained in the subsystems A and B:

    𝒜_A ≡{ĝ∈𝒜| ĝ =ĝ_A⊗1̂}
       𝒜_B ≡{ĝ∈𝒜| ĝ = 1̂⊗ĝ_B }.

The number of elements in these groups is the total number of possible modifications of the state induced by the stars. We can then define the order of these groups as:

    d=|𝒜| = 2^Σ  Σ=number of stars in the systems
        d_A=|𝒜_A| = 2^Σ_A  Σ_A=number of stars fully in A
        d_B=|𝒜_B| = 2^Σ_B  Σ_B=number of stars fully in B.

The presence of stars not fully contained neither in A nor in B makes Σ≠Σ_A +Σ_B. This can be understood from Fig. <ref> (right). The star in the figure depends on a spin in A and on three spins in B. If we perform a trace on one of the two subsystems, its value will then not be completely fixed. To be more precise, the difference in these two quantities is proportional to the length of the boundary between two sub–systems. For the same reason, the deformation operator ĝ can be factorized as:

    ĝ=ĝ_A ⊗ĝ_A/B⊗ĝ_B.

It is useful to introduce the short-hand notation  

    G = ∑_ĝ∈𝒜ĝ∏_i^A/Bσ̂^x_i,

which relates to the operators acting on the A subsystems outside the boundary. 

The sum over the possible modifications inside B can be performed trivially and it gives a term proportional to the order of the group

    ∑_g_B  g_B^'δ( g_B-g_B^') = ∑_g_B1̂= d_B.

The summation over the boundary states, labelled by l, is a little bit more involved but, in general, the states with the same boundary configuration over B are traced out:

    ∑_l,k a_l a_k^∗δ_n_l^B^n_k^B= ( ∑_l_B| a^n_l^B|^2 ) ∑_l_A, k_A a^n_l^A a^n_k^Ae^i(φ_l^A - φ_k^A).

We can then define

    ∑_l_B| a^n_l^B|^2= f(B).

The function f(B) depends just on the boundary of B and on the links between A and B  and it is equal to 1 when it contains all the boundary as far as ∑_l| a^l_n|^2 = 1 for the whole system. From Eq. (<ref>), it is clear that f(B)≤ 1 and it is an increasing function of the argument.


For our future calculations it is useful to work out G|0⟩⟨0|GG̃|0⟩⟨0|G̃:

    ⟨0|GG̃|0⟩ = ⟨0|∑_ĝ∈𝒜ĝ∏_i^A/Bσ̂^x_i ∑_ĝ∈𝒜ĝ'∏_j^A/B'σ̂^x_j |0⟩ .

In order to have non–zero matrix element, we should have

    ∑_ĝ∈𝒜ĝ∏_i^A/Bσ̂^x_i ∑_ĝ∈𝒜ĝ'∏_j^A/B'σ̂^x_j = 1̂.

This condition can be satisfied in two ways:


  * ∏_i^A/Bσ̂^x_i ∏_j^A/B'σ̂^x_j = 1̂. In this case two sets of σ̂^x are identical. Then we should also have ∑_ĝ∈𝒜ĝ∑_ĝ' ∈𝒜ĝ' = 1̂ and there are |d_A| such possibilities.

  * ∏_i^A/Bσ̂^x_i ∏_j^A/B'σ̂^x_j  = Ŝ, where Ŝ = ∏_i^ A/Bσ̂^x_i taken along all the boundary between A and B. In this case, we say that two sets of σ̂^x_i are complementary and together cover all the boundary. Group elements ĝ should then satisfy the next condition:  ∑_ĝ∈𝒜ĝ∑_ĝ' ∈𝒜ĝ' =Ŝ. There are also |d_A| such possibilities.

It is clear, that in any other case (when we have some spins excited along the A/B boundary), this is not the element of group 𝒜 and therefore can not be equal to ∑_ĝ∈𝒜ĝ∑_ĝ∈𝒜ĝ' which leads to zero matrix element. Now we can see that 

    G|0⟩⟨0|GG̃|0⟩⟨0|G̃ = |d_A||0⟩⟨0|G(1̂ + Ŝ).

Now our reduced density matrix can be written in a more elegant way:

    ρ̂_A = d_B f( B )/𝒜×
       ×∑_l,ka^l+k e^i(φ_l-φ_k)_from AĜ∏_i^ℒ_𝒜( 𝒲_i )^e_i|0⟩_A⟨0|_A∏_i^ℒ_𝒜( 𝒲_i )^f_iĜ.

It is also easily seen that
(1̂ + Ŝ)^2 = 2(1̂ + Ŝ).
The expression of the density matrix alone is not enough to compute the entanglement entropy. In order to compute the logarithm of this operator we need to compute all its integer powers.


Starting with the square of ρ̂_A we have:

    ρ̂^2_A   = d_A d_B^2f(B)^2|𝒜|^2∑_k a^2k_=f(A)∑_l,ka^l+k e^iφ_l - iφ_k
       Ĝ[ ∏_i^ℒ_𝒜( 𝒲_i )^e_i|0⟩_A ⟨0|_A ∏_i^ℒ_𝒜( 𝒲_i )^f_i]Ĝ(1̂ + Ŝ) 
       = d_A d_B^2f(B)^2f(A)|𝒜|^2∑_l,ka^l+k e^iφ_l - iφ_k
       Ĝ[ ∏_i^ℒ_𝒜( 𝒲_i )^e_i|0⟩_A ⟨0|_A ∏_i^ℒ_𝒜( 𝒲_i )^f_i]Ĝ(1̂ + Ŝ) 
       =f(B)f(A)d_A d_Bdρ̂_A(1̂ + Ŝ).

We can then proceed to higher powers:

    ρ̂^3_A    = f(B)f(A)d_A d_Bdρ̂^2_A(1̂ + Ŝ)= 
       = 2 ( f(B)f(A)d_A d_Bd)^2ρ̂^2_A.

For the general power of ρ̂_A:

    ρ̂^n_A = ( 2f(B)f(A)d_A d_Bd)^n-2ρ̂^2_A.

Expanding the logarithm in power series it is possible to compute the entanglement entropy as

    S_A = ln( d/2d_A d_B1/f(A)f(B)).

Let us take a look at the expression for the entanglement entropy when the subsystem A is entirely in the bulk. It can be easily done as far as this is just a particular case of the calculation above where the whole boundary is traced out. This would lead to  ∑_l| a^l_n|^2 = 1 already on the stage of calculating the first power of ρ̂_A due to the orthogonality. In the end the expression will be the same but with f(A)f(B) = 1
    S̃_A = ln( d/2 d_A d_B).

Taking into account that d = 2^Σ, d_A = 2^Σ_A, d_B = 2^Σ_B and Σ = Σ_A + Σ_B + Σ_A/B, 
    S̃_A= (ℒ_A/B - 1) ln 2,

which is exactly the result for KTC with PBC.
 Thy we have:

    S_A = S̃_A + ln1 f(A)f(B).

The factor f(A)f(B) is always smaller or equal to  1. It depends just on the portions of boundary present in the subsystems A and B. This contribution leads to the additional term -ln f(A)f(B) which increases the entropy.


 §.§ Entanglement entropy in the bulk

We can now briefly describe the specific case of A completely contained inside the bulk of the system. There, the boundary is completely traced out by the first trace procedure, which then ends up in the reduced density matrix:

    ρ_A = d_B/dĜ|0⟩⟨0|Ĝ,

which gives us the entanglement entropy
S_A =S̃_A
that is exactly the value of the entanglement entropy with PBC. As we explained in the previous sections, the Hilbert spaces of the bulk and the boundary are decoupled. If the system A is completely contained in the bulk, the information about the boundary is then destroyed once the trace over the subsystem B is performed. The density matrix of A is then the same as in the case of PBC, together with the entanglement entropy. The result is as expected as the boundary is entangled.

As a conclusive remark, it is important to notice a significant property of the previous expressions: Eq. (<ref>) and (<ref>) are both completely independent of the value of the couplings. Their values are fixed by the particular geometry we choose for the subsystems A and B.



§ EXCITATION OF THE BOUNDARY THEORY

Hamiltonian (<ref>) has different excitation energies in the bulk and on the boundary due to the different operators acting on these two distinct regions. If we act on the ground state with a generic spin flip, via a σ_i^x operator, it will cost an energy equal to 2( J_m +J_e ) because every spin is shared between two star and two plaquette operators. On the boundary, we have fewer operators acting on our spins, so if we flip the spin at the boundary we need to pay just J_m because just one plaquette operator is acting on this degree of freedom.

Due to the particular structure of the ground state on the boundary, described by 𝒲_i operators, one can see that only two possible excited states are possible: σ_i^x 𝒲_i and σ_i^x σ_i+1^x𝒲_i due to the spin-1/2 algebra. All these excited states are fixed to the points where they are created due to the absence of star operators on the boundary, which allows us to conclude that their dispersion relation is flat.

In order to give our boundary theory non-trivial dynamics, we need to have some dispersion relation. This has been done in different topological error correcting codes as the Wen plaquette and the KTC with different OBC by adding a magnetic field <cit.>.
In our case, we also have a magnetic field along x or y–direction which allows us to flip the spins.

    Ĥ_I = h_x∑_i σ̂^x_i

Using ℋ_I as a perturbation of (<ref>), it is possible to compute the tunneling processes induced by the perturbation, see Ref. <cit.>.
We will now compute the tunnelling between the sites n=i and m=i+R induced by the magnetic field, as represented in Fig. <ref>


Following <cit.>, we can compute the transition |m⟩ =σ_m^x |GS⟩→|n⟩=σ_n^x|GS⟩ using the Gell-Mann-Low theorem <cit.>.

The energy shift is defined as

    E = ⟨n|ĤÛ_I(0,∞)|m⟩ = E_0 + δ E 
       δ E= ⟨n|Ĥ_IÛ_I(0,∞)|m⟩ = ⟨n|Ĥ_I∑_j(Ĥ_IE - Ĥ_0)^j|m⟩

In order to have a well–defined power series expansion, we will assume J_m ≫| h_x|, in this way a perturbative expansion of the previous exact expression is motivated. This is particularly important in order to understand the physics behind these processes. When no magnetic field is activated, the part of the 𝒲_i operator that can be modified by the stars, namely the spin–flip is labelled by α in Eq. (<ref>) and Fig. <ref>. When the magnetic field is activated, every spin flip costs energy. The paths now have different energies depending on their length, this means that we reduce some degeneracy and the shortest paths connecting two points are energetically convenient.

We can start with the nearest–neighbour hopping |i⟩→|i+1⟩ which is clearly i–independent because of the translational invariance of the boundary. These states can be represented using the 𝒲_i ( R ) operators defined before, since they are also the shortest paths connecting two points on the boundary. We then want to compute the transition probability of the process:

    |R=0⟩ = 𝒲_i |0⟩  →  |R=1⟩ = 𝒲_i ( 1 ) |0⟩.

We can compute the energy shift at the leading order as:

    δ E = ⟨R=1|Ĥ_I ( Ĥ_I/E-Ĥ_0)^2 |R=0⟩ = h_x^3/( -2J_m )^2.

An easy way to compute these transitions is using a graphical approach, we draw the lines corresponding to the generic states |R⟩ and for every application of Ĥ_I we flip one spin

    δ E    = ⟨0|𝒲_i𝒲_i+1Ĥ_I(Ĥ_IE - Ĥ_0)^2𝒲_i|0⟩
        = h_x(-2J_m)⟨0|𝒲_i𝒲_i+1Ĥ_I(Ĥ_IE - Ĥ_0)σ̂_x^ασ̂_x^i|0⟩
        = (h_x)^2(-2J_m)^2⟨0|𝒲_i𝒲_i+1Ĥ_Iσ̂_x^βσ̂_x^i𝒲_i|0⟩= (h_x)^3(-2J_m)^2.





Using the same technique, it is possible to evaluate the energy shift for more "long-range" transitions like |R=0⟩→|R⟩ of which the energy shift at the leading order is 

    δ E_R = (h_x)^R+2/(-2J_m)^R+1.




§ DISPERSION RELATION OF THE BOUNDARY THEORY

We want to study the dispersion relation of the boundary theory. Excitations of this theory are composed by σ_i^x operators. We can then use as an effective Hamiltonian the following spin Hamiltonian

    Ĥ_eff = ∑_i,j J_i,jσ̂^x_i σ̂^x_j.

As in Ref. <cit.>, J_i,j is the probability of transitions connecting two spins on the boundary at a distance R=| i -j | given in Eq. (<ref>). The physics of long-range coupling in the Ising model can lead to extremely interesting effects both in equilibrium <cit.> and out-of-equilibrium <cit.> physics. The behavior of the KTC and the stability of topological order in the presence of a magnetic field was already studied in Ref. <cit.>. There, the authors derived an effective description for the bulk Hamiltonian under the perturbation of generic magnetic field. In our set up the magnetic field is chosen to be small enough to do not perturb the bulk but just the boundary degrees of freedom.

We can use the translational invariance of the system to write the dispersion relation of the excitations in Fourier space. We transform our Hamiltonian in Fourier space

    σ_i^x=1/√(ℒ)∑_k e^-  k iσ_k^x.

We can then plug in this expansion into Eq. (<ref>) to get:

    ℋ_eff=1/ℒ∑_i R∑_k q J_R e^(k+q)ie^ qRσ_k^xσ_q^x,

where ℒ is the length of the boundary of our two–dimensional system. The summation over i can be performed and it gives a delta function which enforces the momentum conservation

    ℋ_eff=∑_R ∑_k J_Re^ k Rσ_k^x σ_-k^x.

The previous relation can be used to identify the dispersion relation:

    ϵ_k=∑_R J_R e^ R k.

The realness of the dispersion relation is ensured by the fact that J_R is symmetric under R→ -R

    ∑_R J_R e^ k R = 2∑_R>0 J_R e^ k.

We can then use the Eq. (<ref>)

    J_R=(h_x)^R+2(-2J_m)^R+1=h_x/(-2J_m)^2(h_x/-2J_m)^R.


We can write the dispersion relation using the previous expression and the fact that the sum over R runs from R=1 to R=ℒ/2 which is the maximal distance on the boundary.

The summation can be performed analytically using the λ=h_x/(-2J_m) as a short-hand notation for the perturbation parameter. We get:

    ϵ_k = 2h_x λ∑_R=1^R=ℒ/2(e^ kλ)^R=2h_x λλ e^ k-(λ e^ k)^ℒ/2+1/1-λ e^ k,

the dispersion relation can then be written in the form:

    ϵ_k = 2h_x λ^2cos( k )-λ +λ^ℒ/2+1cos( kℒ/2)-λ^ℒ/2cos( k(ℒ+2)/2)/1+λ^2-2λcos( k).

Taking ℒ sufficiently large and the parameter λ sufficiently small, we can evaluate the previous equation in the thermodynamic limit with respect to the border:

    ϵ_k = 2h_x λ^2 cos(k)-λ/1+λ^2-2λcos(k).
 
An important approximation can be done in the last expression, if λ is sufficiently small, the dominant hopping is the nearest–neighbour hopping, in this case we recover a cosine band

    ϵ_k ≈ 2h_x λ^2 cos(k),

for values of λ not negligible. It is possible to see deviations, meaning that the couplings between more distant sites become important, see Fig. <ref>




§ QUANTUM COMPUTATION WITH THE KTC WITH OPEN BOUNDARY CONDITIONS

In the end, we also want to discuss the question of whether is it possible to use  KTC with OBC for quantum computation. This question is important to discuss as far as this is of the use for the original KTC with periodic boundary conditions and one should understand if and how does the situation changes when one moves to open boundaries. We can start from the considerations on the standard KTC with PBC, where the information is stored in two topologically protected ground states, winding along the two non–trivial paths of the torus. The memory of such a system is protected by an energy gap. The only error that cannot be detected is the one that does not increase energy, meaning that it has to convert one ground state to another one. Since the ground states are paths along the two main directions of the torus, this requires an exponentially large number of ordered spin-flips. Assuming a constant probability to create an error in the system with one spin–flip, the probability to corrupt the information is exponentially small in the size of the system.
 
In our system, the degrees of freedom of the ground states are located at the boundary, given by the 𝒲_i operators. Trivially, one could think about using such operators to store the information one would like to protect. Even though one can naively think that the high complexity and degeneracy of the ground state can make the memory more stable to errors, thus it appears not to be the case. The fundamental reason is that such states are not protected by the topology any more. In fact, to corrupt the information written in our memory, it is sufficient to modify just 3 spins that correspond to creation or destruction of a 𝒲_i operator. This probability is constant and it does not decrease with the system size as in the PBC case. 

Even though, as was seen in the previous section, with the presence of external magnetic field, the edge states are no longer gapless, this still does not make the memory robust. The degenerate ground state splits (similar to Zeeman effect) and acquire an energy gap. One should remember though that the gap is still small due to the fact that it is proportional to magnetic field which is weak. Generally, the gap decreases exponentially with the number of spins excited on the edge. This means that the number of the edge states we can use as q–bits and hence the capacitance of our memory is restricted from the requirement of the gap being sufficiently big. Or, in other words, the level of protection we want to achieve restricts the capacity of the memory we can have.  This is, however, as was said before has nothing to do with the topological protection that the memory has in the KTC with PBC. 

A non-trivial topology can anyway be induced in the system using holes <cit.> and defects <cit.>. Information can then be stored in spin-flips around these points being topologically protected. Planar codes can also be used to do proper quantum computation <cit.> instead of simple storage of information.
Even if our model cannot be directly used as a quantum memory, its utility is clear from the mathematical structure of Eq. (<ref>). This expression can be used to account analytically for the effect of boundary states on the physical quantities and their response to external perturbations.



§ CONCLUSIONS


In this paper, we studied the Kitaev Toric Code with particular open boundary conditions. These conditions restrict any exchange of the excitations between the bulk and the boundary, which are then completely decoupled. We wrote down explicitly the expression for the ground state highlighting the separation between the degrees of freedom in the bulk and the ones localised on the boundary (see Eq. (<ref>)). The presence of the edge modes in open–boundary systems is expected, but our expressions reveal their mathematical form which is not always known. Moreover, the degrees of freedom are located on the boundary but the states themselves traverse deeply into the bulk of the system because of the different operators content of the bulk and the boundary.

We demonstrated that the spatial distribution of the degrees of freedom influences important quantities such as the bipartite entanglement entropy. In fact, the difference between open and periodic boundary conditions can be seen only if the bipartition cuts the boundary. There, the presence of these degrees of freedom gives a positive contribution to the leading order which has an area law behaviour.

Finally, we studied the response of our system to an external perturbation, namely a weak magnetic field. In this situation, the degrees of freedom at the boundary acquire a dispersion relation and they can be described by a generalised Ising model. Using perturbation theory, we calculated the dispersion relation of the model in this limit. Due to the perturbation, the edge modes cannot enter deep into the bulk, but they are mainly confined to the region close to the boundary and consequently take the shortest path between two points.




§ ACKNOWLEDGEMENTS

We thank Kai P. Schmidt and Roman Orus for very helpful discussion. We thank S. Paeckel and K. Harms for useful comments on our manuscript.







10

Haldane1983a
F. D. M. Haldane,
Nonlinear field theory of large-spin heisenberg antiferromagnets:
  Semiclassically quantized solitons of the one-dimensional easy-axis néel
  state,
Phys. Rev. Lett., 50:1153–1156, Apr 1983.

Haldane1988
F. D. M. Haldane,
Model for a quantum hall effect without landau levels:
  Condensed-matter realization of the "parity anomaly",
Phys. Rev. Lett., 61:2015–2018, Oct 1988.

Haldane1983b
F. D. M. Haldane,
Continuum dynamics of the 1-d heisenberg antiferromagnet:
  Identification with the o(3) nonlinear sigma model,
Physics Letters A, 93(9):464 – 468, 1983.

Kosterlitz1974
J. M. Kosterlitz,
The critical properties of the two-dimensional xy model,
Journal of Physics C: Solid State Physics, 7(6):1046, 1974.

Kosterlitz1972
J. M. Kosterlitz and D. J. Thouless,
Long range order and metastability in two dimensional solids and
  superfluids. (application of dislocation theory),
Journal of Physics C: Solid State Physics, 5(11):L124, 1972.

Kosterlitz1973
K. M. Kosterlitz and D. J. Thouless,
Ordering, metastability and phase transitions in two-dimensional
  systems,
Journal of Physics C: Solid State Physics, 6(7):1181, 1973.

Laughlin1981
R. B. Laughlin,
Quantized hall conductivity in two dimensions,
Phys. Rev. B, 23:5632–5633, May 1981.

Kane2005a
C. L. Kane and E. J. Mele,
Quantum spin hall effect in graphene,
Phys. Rev. Lett., 95:226801, Nov 2005.

Kane2005b
C. L. Kane and E. J. Mele.
Z_2 Topological Order and the Quantum Spin Hall Effect,
Phys. Rev. Lett., 95:146802, Sep 2005.

Bernevig2006
A. Bernevig and Shou-Cheng Zhang,
Quantum spin hall effect,
Phys. Rev. Lett., 96:106802, Mar 2006.

Murakami2006
Shuichi Murakami,
Quantum spin hall effect and enhanced magnetic response by spin-orbit
  coupling,
Phys. Rev. Lett., 97:236805, Dec 2006.

Hasan2010
M. Z. Hasan and C. L. Kane,
Colloquium: Topological insulators,
Rev. Mod. Phys., 82:3045–3067, Nov 2010.

Laughlin1983
R. B. Laughlin,
Anomalous quantum hall effect: An incompressible quantum fluid with
  fractionally charged excitations,
Phys. Rev. Lett., 50:1395–1398, May 1983.

Chang2013
Cui-Zu Chang, Jinsong Zhang, Xiao Feng, Jie Shen, Zuocheng Zhang, Minghua Guo,
  Kang Li, Yunbo Ou, Pang Wei, Li-Li Wang, Zhong-Qing Ji, Yang Feng, Shuaihua
  Ji, Xi Chen, Jinfeng Jia, Xi Dai, Zhong Fang, Shou-Cheng Zhang, Ke He, Yayu
  Wang, Li Lu, Xu-Cun Ma, and Qi-Kun Xue,
Experimental observation of the quantum anomalous hall effect in a
  magnetic topological insulator,
Science, 340(6129):167–170, 2013.

Tsui1999
Daniel C. Tsui,
Nobel lecture: Interplay of disorder and interaction in
  two-dimensional electron gas in intense magnetic fields,
Rev. Mod. Phys., 71:891–895, Jul 1999.

Zhang2005
Yuanbo Zhang, Yan-Wen Tan, Horst L. Stormer, and Philip Kim,
Experimental observation of the quantum hall effect and berry's phase
  in graphene,
Nature, 438:201 EP –, 11 2005.

Wen2003a
X-G Wen,
Quantum order from string-net condensations and the origin of light
  and massless fermions,
Phys. Rev. D, 68:065003, Sep 2003.

Wen2003b
X-G Wen,
Quantum orders in an exact soluble model,
Phys. Rev. Lett., 90:016803, Jan 2003.

Kitaev2006
A. Yu. Kitaev,
Anyons in an exactly solved model and beyond,
Annals of Physics, 321(1):2 – 111, 2006.
January Special Issue.

Kitaev2003
A.Yu. Kitaev,
Fault-tolerant quantum computation by anyons,
Annals of Physics, 303(1):2 – 30, 2003.

bravyi1998quantum
S. Bravyi and A. Yu. Kitaev,
Quantum codes on a lattice with boundary,
arXiv preprint quant-ph/9811052, 1998.

Corcoles2015
A. D. Córcoles, Easwar Magesan, Srikanth J. Srinivasan, Andrew W. Cross,
  M. Steffen, Jay M. Gambetta, and Jerry M. Chow,
Demonstration of a quantum error detection code using a square
  lattice of four superconducting qubits,
Nature Communications, 6:6979 EP –, 04 2015.

Beigi2011
Salman Beigi, Peter W. Shor, and Daniel Whalen,
The quantum double model with boundary: Condensations and symmetries,
Communications in Mathematical Physics, 306(3):663–694, Sep
  2011.

Levin2013
Michael Levin,
Protected edge modes without symmetry,
Phys. Rev. X, 3:021009, May 2013.

Yu2013
Jing Yu, Xing-Hai Zhang, and Su-Peng Kou,
Majorana edge states for 𝕫_2 topological orders of the
  wen plaquette and toric code models,
Phys. Rev. B, 87:184402, May 2013.

Legget2013
A. J. Leggett,
The Kitaev models,
Illinnois.edu, 2013.

Hamma2005
P. Zanardi A. Hamma, R. Ionicioiu,
Ground state entanglement and geometric entropy in the kitaev model,
Physics Letters A, 337(1):22 – 28, 2005.

Amico2008
Luigi Amico, Rosario Fazio, Andreas Osterloh, and Vlatko Vedral,
Entanglement in many-body systems,
Rev. Mod. Phys., 80:517–576, May 2008.

Horodecki2009
Ryszard Horodecki, Paweł Horodecki, Michał Horodecki, and Karol
  Horodecki,
Quantum entanglement,
Rev. Mod. Phys., 81:865–942, Jun 2009.

Hamma2005b
A. Hamma, R. Ionicioiu, and P. Zanardi,
Bipartite entanglement and entropic boundary law in lattice spin
  systems,
Phys. Rev. A, 71:022315, Feb 2005.

Srednicki1993
M. Srednicki,
Entropy and area,
Phys. Rev. Lett., 71:666–669, Aug 1993.

Eisert2010
J Eisert, M Cramer, and M B Plenio,
Colloquium,
Rev. Mod. Phys., 82:277–306, Feb 2010.

Kitaev2006topo
A. Yu. Kitaev and Preskill. J.,
Topological entanglement entropy,
Phys. Rev. Lett., 96:110404, Mar 2006.

Levin2006
Michael Levin and Xiao-Gang Wen,
Detecting topological order in a ground state wave function,
Phys. Rev. Lett., 96:110405, Mar 2006.

shi2018characterizing
Bowen Shi and Yuan-Ming Lu,
Characterizing topological orders by the information convex,
arXiv preprint arXiv:1801.01519, 2018.

Kou2009
S-P Kou,
Topological quantum computation via the quantum tunneling effect,
Phys. Rev. Lett., 102:120402, Mar 2009.

Kou2009b
S-P Kou,
Realization of topological quantum computation with planar codes,
Phys. Rev. A, 80:052317, Nov 2009.

Yu2009
J Yu and S-P Kou,
Macroscopic quantum tunneling effect of Z_2 topological order,
Phys. Rev. B, 80:075107, Aug 2009.

GellMann1951
M. Gell-Mann and F. Low,
Bound states in quantum field theory,
Phys. Rev., 84:350–354, Oct 1951.

Belgio2018
Laurens Vanderstraeten, Maarten Van Damme, Hans Peter Büchler, and Frank
  Verstraete,
Quasiparticles in quantum spin chains with long-range interactions,
arXiv preprint arXiv:1801.00769, 2018.

tagliacozzo2012
Thomas Koffel, M. Lewenstein, and Luca Tagliacozzo,
Entanglement entropy for the long-range ising chain in a transverse
  field,
Phys. Rev. Lett., 109:267203, Dec 2012.

Cevolani2015
Lorenzo Cevolani, Giuseppe Carleo, and Laurent Sanchez-Palencia,
Protected quasilocality in quantum systems with long-range
  interactions,
Phys. Rev. A, 92:041603, Oct 2015.

Cevolani2016
Lorenzo Cevolani, Giuseppe Carleo, and Laurent Sanchez-Palencia,
Spreading of correlations in exactly solvable quantum models with
  long-range interactions in arbitrary dimensions,
New Journal of Physics, 18(9):093002, 2016.

Cevolani2017
L. Cevolani, J. Despres, G. Carleo, L. Tagliacozzo, and
  L. Sanchez-Palencia,
Universal Scaling Laws for Correlation Spreading in Quantum Systems
  with Short- and Long-Range Interactions,
ArXiv e-prints, June 2017.

wootton2011bringing
James R Wootton and Jiannis K Pachos,
Bringing order through disorder: Localization of errors in
  topological quantum memories,
Physical review letters, 107(3):030503, 2011.

trebst2007breakdown
Simon Trebst, Philipp Werner, Matthias Troyer, Kirill Shtengel, and Chetan
  Nayak,
Breakdown of a topological phase: Quantum phase transition in a loop
  gas model with tension,
Physical review letters, 98(7):070602, 2007.

Brown2017
Benjamin J. Brown, Katharina Laubscher, Markus S. Kesselring, and James R.
  Wootton,
Poking holes and cutting corners to achieve clifford gates with the
  surface code,
Phys. Rev. X, 7:021029, May 2017.

Brown2013
Benjamin J Brown, Stephen D Bartlett, Andrew C Doherty, and Sean D Barrett,
Topological entanglement entropy with a twist,
Physical review letters, 111(22):220402, 2013.

bombin2007optimal
H Bombin and Miguel A Martin-Delgado,
Optimal resources for topological two-dimensional stabilizer codes:
  Comparative study,
Physical Review A, 76(1):012305, 2007.

bombin2007homological
Hector Bombin and Miguel A Martin-Delgado,
Homological error correction: Classical and quantum codes,
Journal of mathematical physics, 48(5):052105, 2007.

nigg2014quantum
Daniel Nigg, Markus Mueller, Esteban A Martinez, Philipp Schindler, Markus
  Hennrich, Thomas Monz, Miguel A Martin-Delgado, and Rainer Blatt,
Quantum computations on a topologically encoded qubit,
Science, 345(6194):302–305, 2014.





