HOWTO /
ToricCodeIn this tutorial, we will introduce the toric code model, walk through the analytic solution of the model, and discuss performing calculations of this model using the tookit. IntroductionThe toric code model was introduced by Kitaev in 2003 as a proof of concept to show that physical materials with topological order could be used store quantum information in a way which is robust to errors [1]. This is due to the fact that topologically ordered systems have a ground state degeneracy which is protected from breaking under arbitrary perturbations to the Hamiltonian which do not cause a phase transition. This ground state degeneracy is directly related to topology of the system, and on a torus, there are four ground states: this corresponds to two computational qubits. As we shall see, the ground state degeneracy is also related to the elementary excitations of the system, which are anyons (that is, they do not behave as normal bosons or fermions). By creating a pair of anyons on the ground state, and moving them around a noncontractible loop on the torus and annihilating them, we can perform operations on the computational qubits. Hence, the ground states cannot be transformed into each other using only local operations (as moving excitations around a noncontractible loop is a global operation), and so the computational qubits are protected from errors by the underlying physical properties of the system if it is sufficiently cooled. However, the toric code model is not very physical (as it has fourspin interactions) and Abelian topological order like in the toric code is considered less useful for quantum computation than the nonAbelian case, since universal quantum computation can be realised by braiding nonAbelian anyons around each other, whereas moving the Abelian anyons around the torus cannot be used to realise every quantum gate. Nevertheless, because it is (relatively) intuitive to understand, and is exactly solvable, the toric code model is still quite useful as an introduction to the ideas of topological order and anyons. The toric code modelThe toric code model is defined on a square lattice, where a spin1/2 particle is placed on each edge of the lattice (not on the vertices), see Fig. 1. For each vertex {$s$} and face {$p$} in the lattice, we define the following star and plaquette operators {$$A_s = \prod_{i\in s} \sigma^x_i, \qquad B_p = \prod_{i\in p} \sigma^z_i,$$} where {$\sigma^\alpha_i$} is the Pauli {$\alpha$}matrix acting on site {$i$} and the products run over the four spins adjacent to that vertex or face. The toric code Hamiltonian is then given by {$$H = J_e \sum_s A_s  J_m \sum_p B_p,$$} where {$J_e$} and {$J_m$} define the energy scale of the system. Since each {$A_s$} and {$B_p$} operator commutes with every other {$A_s$} and {$B_p$} operator (you can prove this yourself), they can all be simultaneously diagonalised. To minimise the energy, the state must be an eigenstate of each {$A_s$} and {$B_p$} with eigenvalue {$+1$}, and so the ground states will be in the subspace defined by the property {$A_s = B_p = +1$}. ExcitationsIn an excited state, we must then have that {$A_s = 1$} or {$B_p = 1$} for at least one star or plaquette: hence, we identify a star with {$A_s = 1$} as an elementary excitation with energy {$2J_e$}, which we call an electric charge or an {$e$} particle, and a plaquette with {$B_p = 1$} as a magnetic vortex or {$m$} particle with energy {$2J_m$} (the names come from an analogy with a {$Z_2$} lattice gauge theory). We can create a pair of {$e$} or {$m$} particles in the lattice by applying a {$\sigma^z$} or {$\sigma^x$} operator, respectively, onto the ground state, as shown in Fig. 2. Each particle is its own antiparticle, so applying a Pauli matrix on a spin adjacent to an excitation will destroy that excitation. Hence, we can create two particles at a distance by applying 'strings' of Pauli operators to the ground state {$$Z[t] = \prod_{i\in t} \sigma^z_i, \qquad X[t'] = \prod_{i\in t'} \sigma^x_i,$$} where {$t$} ({$t'$}) is a string of sites with adjacent vertices (faces). This can be physically interpreted as creating a pair of particles from the vacuum and then moving one along the string: see Fig. 3. Similarly, applying a closed loop of Pauli matrices can be interpreted as creating a pair of particles from the vacuum, moving one particle along the loop, and then annihilating them again. For a contractible loop of {$\sigma^z$} ({$\sigma^x$}) operators, this is equal to the product of all plaquette (star) operators inside of the loop, and so the expectation value of this operation will be the parity of {$m$} ({$e$}) particles inside of the loop (in Fig. 3, the expectation value of {$Z[l]$} will be {$1$}). From this, we can see that these excitations have anyonic particle statistics: for ordinary bosons and fermions, moving one particle around another will not affect the phase of the state, but with the Abelian anyons in the toric code, if we move an {$e$} particle around an {$m$} particle, we will gain a phase of {$1$}. (We note that while the {$e$} and {$m$} particles have bosonic selfstatistics, i.e. they do not change phase upon exchanging the position of two {$e$} or two {$m$} particles, they are anyons because their mutual statistics is nontrivial. Furthermore, if we fuse and {$e$} and {$m$} particle, we obtain an {$\epsilon$} particle, which has fermionic selfstatistics, which is not what normally happens when we combine two bosons!) We note that you cannot get a single isolated excitation with only local operations: the number of each type of particle will always be even. The only way to get a single excitation when working in the thermodynamic limit is to create a pair and move one of the particles 'off to infinity'. Ground state degeneracyNow that we understand the excitations of the toric code model, we will discuss the topologicallydependent ground state degeneracy. Suppose we write a ground state {$\Psi\rangle$} of the toric code in terms of the basis {$\{\psi_i\rangle\}_i$} formed by tensor products of the eigenstates of the {$\sigma^z$} operators {$\uparrow\rangle$} and {$\downarrow\rangle$}, then any ground state must be a superposition of 'vortexfree' basis states where there are an even number of up and down spins around each plaquette (i.e. {$B_p\psi_i\rangle = \psi_i\rangle$} for each {$p$}), so {$$\Psi\rangle = \sum_i c_i \psi_i\rangle, \quad \text{where} \quad B_p\psi_i\rangle = \psi_i\rangle \quad \forall p.$$} Now if {$\psi_j\rangle = A_s \psi_i\rangle$} for some star {$s$}, for {$\Psi\rangle$} to be a ground state, we need {$A_s \Psi\rangle = \Psi\rangle$}, and so we need that {$c_i = c_j$}. Now, for any two vortexfree states {$\psi_i\rangle$} and {$\psi_j\rangle$}, we can get from {$\psi_i\rangle$} to {$\psi_j\rangle$} by applying some number of closed loops of {$\sigma^x$} operators {$X[l']$} (try to convince yourself that this is true). If all of these loops are contractible, then they can be written as the product of {$A_s$} operators acting inside of the loops, and then we must have that the coefficients {$c_i$} and {$c_j$} are equal. If the system is on the surface of a manifold such as a sphere or an infinite plane, all loops are contractible, and so all the coefficients {$c_i$} must be equal in the basis decomposition, and there is therefore only one unique ground state. For a manifold such as a torus or an infinite cylinder, we have loops that are noncontractible, and so applying the operator {$X[l']$} around these loops cannot be expressed in terms of a product of star operators, and the ground state is not necessarily unique. On a torus, there are exactly two independent noncontractible loops, and we can choose {$l_x$} and {$l_y$} in Fig. 4 to be these loops, and so we can label the ground states by the expectation values of {$X[l_x]$} and {$X[l_y]$}. Since {$X[l_x]$} and {$X[l_y]$} cannot be decomposed into star operators, they can be either {$+1$} or {$1$} on the ground state. Therefore, there is a fourfold ground state degeneracy on the torus. To transform from one ground state to the other we can apply the {$Z[l_x]$} and {$Z[l_y]$} operators. This is since, for a contractible loop {$l$}, {$X[l]$} measures the parity of particles inside of the loop by counting the parity of {$Z[t]$} strings crossing the loop {$l$}: similarly, for a noncontractible loop {$l$}, {$X[l]$} changes sign if we act with a string {$Z[t]$} crossing the loop exactly once (or an odd number of times). Hence, {$Z[l_x]$} transforms from the ground state with {$X[l_y] = +1$} to {$1$}. And so, these ground states are physically interpreted as having the presence ({$1$}) or absence ({$+1$}) of an {$m$} particle flux crossing the loops {$l_x$} and {$l_x$}, as shown in Fig. 5. Minimally entangled statesThis particular choice of basis for the ground space is helpful for understanding how the topological ground state degeneracy arises, but there is another basis which is more important to use when simulating topologically ordered systems using MPS techniques. The first thing to note is that when we wish to model a 2D system using MPSs, we place the system on the surface of a cylinder where the circumference of the cylinder {$L_y$} is very small (usually less than {$10\sim 20$} sites), while the length of the cylinder {$L_x$} will be much larger than {$L_y$}, if not infinite. Now if we look at the singular value spectrum at any bond in our MPS, this will correspond to a Schmidt decomposition of our system where we cut the system into halves along the circumference. Because of the nature of DMRG and other MPS ground state search algorithms, ground states where the Schmidt values decay faster will be favoured over states where they decay slower: in other words, ground states with a smaller entanglement entropy are favoured, the socalled minimally entangled states (MESs). (This is since if we truncate to some fixed bond dimension, a minimally entangled ground state will be represented with greater fidelity than a more entangled ground state.) It is then the basis of MESs which are favoured when finding the ground states of a topologically ordered, or even a spontaneoussymmetrybroken system. For example, consider the transversefield Ising model in the symmetrybroken phase: at zero field, there are two minimally entangled ground states: {$\Psi_\uparrow\rangle = \cdots\uparrow\uparrow\uparrow\uparrow\uparrow\cdots\rangle$} and {$\Psi_\downarrow\rangle = \cdots\downarrow\downarrow\downarrow\downarrow\downarrow\cdots\rangle$}, both of which have bond dimension 1 when written as an iMPS. Any superposition {$\alpha\Psi_\uparrow\rangle + \beta\Psi_\downarrow\rangle$} is also a ground state, but will have a bond dimension of 2. Similarly, for a nonzero field below the critical value, we will still have two minimally entangled ground states {$\Psi_\uparrow\rangle$} and {$\Psi_\downarrow\rangle$} where the Schmidt values decay the fastest, and any superposition {$\alpha\Psi_\uparrow\rangle + \beta\Psi_\downarrow\rangle$} will have slower decaying Schmidt values, and will not be approximated as well if we try to represent it as a finite bond dimension iMPS. For the toric code, the minimally entangled ground states are the simultaneous eigenstates of the {$X[l_y]$} and {$Z[l_y]$} operators, which are the states that can be interpreted as containing an anyonic flux going down the length of the cylinder (for justification as to why these are the MESs, see Ref. [2]): these states are shown in Fig. 6. Since each state corresponds to a different anyonic flux going down the cylinder, together, the set of MESs can provide information about many properties of the excitations. Considerations for other topologically ordered modelsMany of the ideas presented above for the toric code extend intuitively to other anyonic models, such as the ground state degeneracy on a torus corresponding to the number of species of anyons. One important caveat is that many other models with possible emergent topological order, such as the Heisenberg model on a frustrated lattice, will not be as numerically tractable as the toric code. While each of the MESs can be found easily for the toric code, for more physically realistic models, it can be a much more difficult task to find all of these states for a specific cylinder width, since the ground state search algorithm may be biased towards some of them more than others. As a result, actually using all of the MESs to investigate the topological order of a system can turn out to be very difficult numerically. Using the toolkitGround statesTo generate a lattice on a width4 cylinder with the toolkit, use the command $ toriccode w 4 o lattice The unit cell of this lattice corresponds to a ring of sites around the circumference, alternating between sites on vertical and horizontal links, as shown in Fig. 7.
As with other models in the toolkit, local operators are specified using an index of the form $ mpidmrgs3e H "lattice:H_starH_plaq" w psi m 16x100 bf (Note that for a cylinder of width {$L_y$}, we only need a maximum bond dimension of {$2^{L_y}$}, 16 in this case, to represent the ground state exactly.) To check which MES this wavefunction is, we can measure the expectation values of the loops around the circumference {$X[l_y]$} and {$Z[l_y]$}, which are implemented as the $ mpexpectation psi "lattice:WX(0)" $ mpexpectation psi "lattice:WZ(0)" In our case, we obtain {$+1$} for both values, although in practice, the obtained state can be in another sector. To transform this state to another MES in a different flux sector, we can apply the {$X[l_x]$} operator, which can be implemented as $ mpiapply "lattice:prod_unit(X(0)[0])" psi psi2 In our case, we can confirm that $ mpiapply "lattice:prod_unit(Z(0)[1])" psi psi3 $ mpiapply "lattice:prod_unit(Z(0)[1])" psi2 psi4 Virtual symmetry operatorsThe different MES sectors can also be distinguished by how symmetry operators, such as translation and reflection, act on the virtual level of the MPS (i.e. on the bonds of the MPS). In particular, the application of a symmetry operator {$U$} on the physical legs of an MPS is equivalent to inserting operators {$\tilde{U}$} in the virtual level of the MPS: these virtual symmetry operators {$\tilde{U}$} form a projective representation of the symmetry group (for a discussion of this phenomenon in the context of symmetryprotected topological order, see Ref. [3] and this article on this site; Ref. [4] discusses this in the context of topological order in the triangular Heisenberg model). For example, we can look at the symmetry operators $ mpauxalgebra w psi l lattice Ry TyPi which we see is {+1} for all states but {$\Psi_\epsilon\rangle$}, which has {$1$} instead. Excitations... References[1] A. Yu. Kitaev, Faulttolerant quantum computation by anyons, Ann. Phys. (NY) 303, 2 (2003), doi:10.1016/S00034916(02)000180, arXiv:quantph/9707021.
[2] Y. Zhang et al., Quasiparticle statistics and braiding from groundstate entanglement, Phys. Rev. B 85, 235151 (2012), doi:10.1103/PhysRevB.85.235151, arXiv:1111.2342.
[3] F. Pollmann and A. M. Turner, Detection of symmetryprotected topological phases in one dimension, Phys. Rev. B 86, 125441 (2012), doi:10.1103/PhysRevB.86.125441, arXiv:1204.0704.
[4] S. N. Saadatmand and I. P. McCulloch, Symmetry fractionalization in the topological phase of the spin{$\frac{1}{2}$} {$J_1$}−{$J_2$} triangular Heisenberg model, Phys. Rev. B 94, 121111(R) (2016), doi:10.1103/PhysRevB.94.121111, arXiv:1606.00334.
