Breakdown of a perturbed Z_N topological phase

We study the robustness of a generalized Kitaev's toric code with Z_N degrees of freedom in the presence of local perturbations. For N=2, this model reduces to the conventional toric code in a uniform magnetic field. A quantitative analysis is performed for the perturbed Z_3 toric code by applying a combination of high-order series expansions and variational techniques. We provide strong evidences for first- and second-order phase transitions between topologically-ordered and polarized phases. Most interestingly, our results also indicate the existence of topological multi-critical points in the phase diagram.


Introduction
The concept of topologically-ordered quantum matter has been introduced by Wen in the context of high-temperature superconductivity [1,2], and is crucial to characterize fractional quantum Hall states and topological insulators. The most striking property of topologically-ordered quantum phases is their dependence on non-local properties of the system. As a consequence, such phases cannot be characterized by a local order parameter, so that the celebrated Landau's symmetry-breaking theory cannot be used.
More recently, topological order has attracted great interest in the field of quantum information due to its weak sensitivity to any local perturbation [3][4][5]. Indeed, non-local degrees of freedom associated with this exotic order have been shown to be (topologically) protected against local sources of decoherence. This key idea is at the heart of topological quantum computation [3,6]. It is thus of importance to quantify precisely this protection when perturbations are added. One prominent example where the effect of additional perturbations has been extensively discussed is the case of Kitaev's toric code [3]. The toric code is an exactly solvable two-dimensional quantum spin model with a Z 2 spin-liquid ground state possessing gapped (Abelian) anyonic excitations. It can be considered as one of the simplest models displaying topological order. Apart from the influence of temperature [7][8][9][10][11] and of disorder [12,13], several works have investigated the effect of an external magnetic field in this model [14][15][16][17][18][19][20]. Interestingly, a very rich phase diagram containing first-and secondorder phase transitions, multi-criticality, self-duality, and dimensional reduction has been found. Similar phase transitions out of topologically-ordered phases have been studied in the context of the Levin-Wen model [21][22][23].
Actually, the toric code model can be defined for any discrete Abelian or non-Abelian group [3,24,25]. In this work, we present an extension of this model to Z N degrees of freedom [26], which reduces to the conventional toric code for N = 2. Although excitations are still Abelian, qualitative and quantitative analyses of this model in the presence of local perturbations provide some insights in the understanding of topological phase transitions when more complex degrees of freedom are involved.
The paper is organized as follows : in section 2, we describe an extension of Kitaev's toric code from Z 2 to Z N degrees of freedom. Our starting point is Wen's plaquette model [27] for which a generalization to Z N can be written down rather easily. The Z N plaquette model is then mapped onto a Z N toric code which is shown to display topological order and anyonic statistics.
To study the robustness of the topological order, we analyze the influence of local perturbations in section 3. For N = 2, such perturbations correspond to a uniform magnetic field. If the perturbations are strong enough, one expects conventional polarized phases that are not topologically ordered. As a consequence, a phase transition between the topological and the polarized phases must occur. Analyzing the breakdown of the topological phase is a very challenging problem for general N . Nevertheless, as we shall see for special kinds of perturbation, either exact mappings onto already known models exist or the model displays self-duality and dimensional reduction.
In section 4, we focus on the case N = 3 and probe the robustness of the Z 3 topological phase for simple perturbations. To this end, we use perturbative continuous unitary transformations (pCUT) and a variational approach based on infinite projected entangled pair states (iPEPS). Note that this combined pCUT+iPEPS method has been already used successfully for the standard toric code (N = 2) in an arbitrary magnetic field [19]. After a brief discussion of both methods and their combination, results are presented for two types of perturbation that can be viewed as natural generalizations of the N = 2 toric code either in a parallel field or in a transverse field. In the simplest case where the perturbation commutes with local charge (or flux) operators, we establish an exact mapping, valid at low-energies, onto a three-state clock model in a transverse field. For a ferromagnetic coupling, this model is known to display a weakly first-order transition (see for instance [28]). However, the perturbation considered here also leads us to study the antiferromagnetic three-state clock model in a transverse field which, to our knowledge, has never been discussed in the literature. In this work, we find evidence for a second-order transition in this system (whose universality class still remains to be accurately determined). We then discuss the counterpart of an arbitrary parallel magnetic field in the N = 2 model. For such a perturbation, we obtain a rich phase diagram containing first-and second-order phase transition lines that form the boundary of the topological phase. Finally, we study the transverse-field problem that diplays dimensional reduction and self-duality, as the N = 2 model [18]. Conclusions and perspectives are drawn in section 5.

Construction of the generalized Z N toric code
In the following, we first introduce the Z N -generalization of the plaquette model introduced by Wen [27]. The main reason to do so is that the plaquette model with Z N degrees of freedom arises rather naturally from the Z 2 conventional counterpart and it is very simple to write down. Afterwards we perform a mapping to a generalized toric code model with Z N -Abelian anyons [26].

Z N plaquette model
We consider N -state degrees of freedom located on the sites i of a square lattice whose unit-cell vectors n 1 and n 2 are shown in figure 1 (left). The associated orthonormal states are denoted by |q i , where q ∈ Z N . Next, let us define the operators Z i and X i as where ω = e 2iπ N . These unitary operators reduce to the conventional Pauli matrices σ z i and σ x i for N = 2. On the same site, both operators obey the important "commutation relation" (Weyl algebra) which generalizes the well-known anticommutation relation σ x i σ z i = −σ z i σ x i of Pauli matrices. They obviously commute when acting on different sites.
The Hamiltonian of the Z N plaquette model is then defined by where W p = Z D X R Z U X L . Sites D, R, U, and L correspond to the four sites (i, i + n 1 , i + n 1 + n 2 , i + n 2 ) of an elementary plaquette p of the square lattice (see figure 1).
Using (1) and (2), it is easy to check that L R U D n 2 n 1 i p for all p and p , so that all plaquette operators commute with the Hamiltonian. Furthermore, it is important to note that these operators obey W N p = 1. In other words, W p 's are Z N conserved quantities whose eigenvalues are simply {ω q , q ∈ Z N }. As for the standard Z 2 plaquette model [27], this property ensures exact solvability of H plaquette . In the following, we map the plaquette model onto a generalized toric code [3,26] as was already done for the Z 2 case [29,30]. The main advantage of the toric code is that the quantum statistics of the elementary excitations are simpler to identify. Additionally, it allows one to adopt an ad hoc language reminiscent of lattice gauge theories with Z N degrees of freedom (see, e.g., [31,32] [3], we introduce a translationally invariant and four-colored lattice as depicted on the right side of figure 1, where the N -state degrees of freedom of the plaquette model are placed on the bonds of a square lattice. We then perform the following local, unitary transformations : As a consequence, the plaquette operators W p become different on the red and blue stars s and on the cyan and pink plaquettes p as illustrated in figure 1. To strengthen the analogy with the conventional toric code, we relabel the various operators as follows : Let us underline that the four-coloring of the lattice is mandatory if one wants to define B p with Z (or Z † ) operators only. Nevertheless, other choices with smaller units cells are possible. Finally, we obtain the Hamiltonian of the Z N toric code Let us note that this model was already introduced in [26] but differs from the Z N toric code discussed in [3] which involves projectors P s and P p instead of A s + A † s and B p + B † p (see next section for definitions). However, since both models are equivalent for N = 2, 3, we shall (abusively) call them toric codes.
2.2.2. Ground states and topological degeneracy A direct consequence of (4) is that all A s and B p operators commute with H TC . Thus, the ground-state energy (per site) e 0 is simply obtained by choosing the (possibly degenerate) minimal eigenvalue of the local operators −J(A s + A † s ) or −J(B p + B † p ), namely, e 0 = −2J cos(2πk/N ). For J > 0 and for any N , the ground state is unique and obtained for k = 0. However, for J < 0, the ground state is unique for N even (in this case, one chooses k = N/2) but it is infinitely-many degenerate for N odd since, locally, one can choose k = (N ±1)/2. In the following, we will only consider the simplest case J > 0.
Nevertheless, there are subtleties since the ground-state degeneracy also depends on the surface's topology as we shall now see on two simple examples. Let us first consider an infinite open plane for which no constraint on A s 's and B p 's exists. In this case, the ground state is unique and can be built as where N is a normalization constant and Operators P s (P p ) project on subspaces with eigenvalue 1 of the corresponding A s (B p ). The reference state |ref can be chosen arbitrarily provided it leads to |gs = 0. For instance, one may choose the fully-polarized state |ref = i |0 i that already fulfills P p |ref = |ref for all plaquettes p.
Next, let us consider the Z N toric code on a torus. In this case, there are two constraints so that the number of independent eigenvalues of A s and B p operators is reduced by two. However, as in the Z 2 toric code, there exist conserved loop operators that can be chosen as : where C k with k ∈ {1, 2, 3, 4} are the non-contractible loops of the torus depicted in figure 2 (left). All these operators commute with H TC , but only two of them can be chosen independently since [Z µ , X µ ] = 0 with µ ∈ {1, 2}. These two additional conserved quantities maintain the exact solvability of H TC . Concretely, if we choose to label states with the eigenvalues z µ = {ω q } (q ∈ Z N ) of the operators Z µ , we find that there exist N 2 ground states |gs, z 1 , z 2 that can be written as where N is a normalization constant. More generally, following [3,26], one can show that for a compact surface with genus g, each eigenstate is N 2g -degenerate (at least) so that the system is indeed topologically ordered.
2  (13) and (14). An example of a counter-clockwise braiding contour for moving a charge qs initially located on a red star s around a flux qp located on a cyan plaquette p is represented in white (see text for explanations). The black circle locates the site (numbered 8) where the crossing between the braiding contour and the string Sp of the flux occurs.

Excitations and statistics
Excitations of the toric code correspond to states which violate the condition that all eigenvalues of the A s or the B p operators are equal to 1. In other words, an elementary particle is a charge q on star s (a flux q on plaquette p) corresponding to an eigenvalue ω q , with q ∈ Z N and q = 0, of A s (B p ) (remember that q = 0 defines the ground state). As A s and B p operators commute with H TC , charges and fluxes are static excitations. Moreover, the form of the Hamiltonian implies that the energy of a many-particle state is simply the sum of the single-particle energies. In other words, charges and fluxes do not interact. In what follows, we give the explicit construction of single-particle states for general N following the detailed construction given in [30] for N = 2. As was already the case in the Z 2 toric code, there is no local operator creating a single excitation. This is illustrated in figure 3 where one can see that an X i (Z i ) operator creates two excitations on neighboring plaquettes (stars).
However, for an infinite plane with open boundary conditions, it is possible to consider single-particle states. Indeed, in such a system, a single excitation can be obtained by first creating a pair of charges (fluxes) and by taking one of the particles to infinity, at least in principle. It is clear that this is rather a gedankenexperiment that cannot be implemented practically (for instance, in a computer). Nevertheless, one may always consider a state where the two particles originating from the elementary pair-creation process are so distant that they are eventually independent. Thus, a one-particle state |q α with charge (flux) q on α = s (p) can be defined, for instance, X Z Figure 3. Illustration of the action of operators X i (left) and Z i (right) on an eigenstate of As and Bp operators. A diamond on a star or a plaquette, with a ± sign denotes a multiplicative change by ω ±1 of the corresponding eigenvalue. The behavior depends on the site's color. as where the semi-infinite strings S α are displayed in figure 2 (right).
We shall now show on a specific but sufficiently general example that charges and fluxes obey mutual anyonic statistics. Let us consider an eigenstate, denoted by |ψ , with a charge q s and a flux q p at positions shown in figure 2 (right). One can braid the charge around the flux along the counter-clockwise oriented white path drawn in this figure, by acting on |ψ with the operator O = Z qs 10 . . . Z qs 7 (Z qs 6 ) † Z qs 5 . . . Z qs 1 . This can be checked from the action of Z i operators shown in figure 3 (right). Furthermore, from the definition of B p operators given in (6), this operator O is nothing but the product of all B † p qs operators for all plaquettes encircled by the braiding contour.
Given that |ψ is an eigenstate of all plaquette operators with eigenvalue 1, except for the plaquette p where the flux is located and for which B p |ψ = ω qp |ψ , one gets O |ψ = ω −qsqp |ψ . This non-trivial braiding phase is the signature of the mutual (Z N ) anyonic statistics between charges and fluxes. It is similar to an Aharonov-Bohm phase, which explains the terminology of charges and fluxes employed to describe the excitations. The same argument allows one to show that a braiding of a charge (flux) around a charge (flux) leads to a trivial phase which is reminiscent from the bosonic statistics of charges (fluxes). In addition, it is clear that there is a hard-core constraint since one cannot create two particles on the same star or on the same plaquette.
Let us end this discussion by showing that the phase can be obtained in another, complementary way. One can write explicitely the state |ψ as Then, one can compute O |ψ by commuting O with the operators appearing in the above expression using (2) and O |gs = |gs since the ground-state is flux-free. For the particular braiding represented in figure 2 (right), the non-trivial phase will appear from the commutation Z qs 8 X qp 8 = ω −qsqp X qp 8 Z qs 8 at site number 8 where the braiding contour and the string of the flux S p intersect.

Z N toric code in the presence of local uniform perturbations
In this section, we first define the general local perturbations for the Z N toric code model. The presence of such perturbations destroys the exact solvability of the toric code model. As mentioned in the introduction, a phase transition has to take place when the perturbation strength increases since, for J = 0, the ground state is fully polarized and thus not topologically ordered. Thereafter, we focus on special examples that allow for a detailed investigation of this transition for general N .

General structure of the perturbation
Here, we shall only consider perturbations which act locally on a site i. A basis of the space of local unitary operators can be conveniently written down in terms of the operators X i and Z i and their powers. Hermitian combinations of these operators lead to the following general form of local uniform perturbations where h l,m ∈ C, and (l, m) ∈ Z 2 N with 0 l N/2. The action of the operator X l i Z m i is displayed in figure 4. As can be inferred from this figure, the perturbation violates the local conservation of charges (fluxes) when m = 0 (l = 0) since it induces creation and annihilation of pairs of excitations as well as hopping processes. However, it violates neither the conservation of the total charge s q s modulo N , nor the conservation of the total flux p q p modulo N [we recall that a star s (a plaquette p) is said to carry a charge q s (a flux q p ) if A s (B p ) has eigenvalue ω qs (ω qp ), with q s and q p being defined modulo N .] Figure 4. Illustration of the action of the operator X l i Z m i on an eigenstate of As and Bp operators. A diamond on a star or a plaquette, with a value k = ±l, ±m, denotes a multiplicative change by ω k of the corresponding eigenvalue. The behavior depends on the site's color.
In fact, these conservation rules can be more conveniently rewritten as conservation rules of "unphysical" total charge and flux, belonging to Z and not to Z N (they will also prove to be useful later on, see section 3.3). Let us denote the "unphysical" charge (flux) at star s (plaquette p) as q s ( q p ). They are related to the true/physical charge (flux) via the following relations q s = q s mod N (q p = q p mod N ). They can be defined non ambiguously by "fixing a gauge" as we explain now (having in mind a purely computational perspective). Given one eigenstate of H TC , it is possible to choose any value for the "unphysical" charges and fluxes (provided they yield the correct physical charges and fluxes). Then, the full Hamiltonian H TC + l,m H l,m can be studied in the subspace of eigenstates of H TC , spanned by the repeated action of the perturbation. Each of these states can be assigned unique values of q s and q p using rules shown in figure 4. It should be clear that s q s and p q p are the same for all states in the generated subspace, which means the Hamiltonian is block diagonal and the total charge s q s and total flux p q p are conserved integers.

Simplest cases
Let us first discuss the case l = 0 for which the perturbation acts only non trivially on charges (which are no longer locally conserved) whereas fluxes remain static gapped excitations that might be seen as sources of Aharonov-Bohm-like phases for the moving charges. Of course, the respective roles of fluxes and charges are exchanged if l = 0 and m = 0.
In this work, we are interested in transitions between the topological phase (existing for small enough perturbations) and the polarized phase expected for J = 0. To get a first idea about the associated physics, let us consider the simplest perturbation corresponding to (l = 0, m = 1), for which the Hamiltonian of the perturbed toric code reads where we set h 0,1 = h Z ∈ R . Since we are interested in the low-energy properties and since [H 0,1 , B p ] = 0 for all p's, we only consider the flux-free subspace where all B p 's have eigenvalue 1, in which the energy of the fluxes is minimal. Then, following the procedure discussed in [14,15] for the special case N = 2, let us denote by |q s the eigenstates of A s 's with eigenvalues ω q . From figure 3 (right), one can check that the action of Z i + Z † i on a site i located between two neighboring stars s and s is equivalent to X s X † s + X † s X s where, following (1), we introduce the operator X s defined by X s |q s = |q − 1 s . Thus, defining Z s = A s (such that Z s |q s = ω q |q s ), one can map Hamiltonian (17) onto the N -state clock model in a transverse field [33] where s, s denotes nearest-neighbor stars s and s . The term −2JN p arises from the replacement of all B p -operators by their eigenvalue 1 (N p denotes the total number of plaquettes). It is important to stress that this mapping preserves neither the degeneracies of the energy levels (hence the topological order) nor the quantum statistics. However, the zero-temperature phase diagrams of H TC + H 0,1 and H clock are exactly the same. Let us remark that the coupling term h Z in (18) stems from the local perturbation in (17) that can be either positive or negative. When h Z > 0, the coupling between stars in H clock is ferromagnetic whereas h Z < 0 leads to antiferromagnetic interactions. This distinction is irrelevant for even N since, in this case (and for a bipartite lattice), one can always perform local unitary transformations which map H clock onto a ferromagnetic model. By contrast, for odd N , one must distinguish between both signs that may lead to various types of transitions.
Unfortunately, few results are available in the literature concerning the twodimensional quantum clock model in a transverse field except for N = 2 (Ising model) where a second-order transition occurs, for N = 3 (Potts model) where a weakly firstorder is expected for h Z > 0 (see for instance [28]), and for N = 4 where the model is equivalent to two decoupled transverse-field Ising model [23]. In such a context, the second-order transition found in section 4 for N = 3 and h Z < 0 opens some interesting perspectives.

Self-duality
The Z 2 toric code in a transverse field is known to be self-dual [18,34]. We shall now show that this property still holds for a general value of N , provided one chooses particular values of l and m. Let us consider the Hamiltonian H TC + H l,m [see (7) and (16)]. This Hamiltonian will be self-dual if its spectrum is symmetric (up to degeneracies) under the exchange J ↔ |h l,m |. Roughly speaking, this will be ensured provided the "roles" of A s and B p operators can be played by the operators X l i Z m i , once stars and plaquettes are exchanged with sites, i.e., when considering the dual lattice illustrated in figure 5.
In the limiting cases, J = 0 and h l,m = 0, self-duality imposes that spectra of H TC and H l,m are the same (up to degeneracies). Setting h l,m = |h l,m |e iφ l,m , this means that e iφ l,m X l i Z m i must have the same spectrum as A s (or B p ). Since this spectrum is {ω q , q ∈ Z N }, this leads to two constraints : whereas the second one can be rephrased as The third and most stringent condition is that the action of H l,m on eigenstates of H TC , i.e., of all A s and B p operators (see figure 4), is the same as the action of H TC on eigenstates of H l,m , i.e., of all X l i Z m i operators (see figure 6). After noticing that the numbers of minus signs, as well as their exact positions, can be made the same in both figures by exchanging the roles of X l i Z m i and its hermitian conjugate on magenta and green sites (or on yellow and brown sites), this last condition reads For N = 2, one recovers the model studied in [18] (up to a factor of 2). Indeed, the only solution to the above equations is l = m = 1 and φ 1,1 = ±π/2. Thus, the perturbation reads ±2h y i σ y i , while the toric code Hamiltonian can be simplified to −2J( s A s + p B p ) since the star and plaquette operators are Hermitian for N = 2. Note that the present analysis also shows that the sign of the field term is irrelevant.
To conclude this discussion, let us show that self-duality in these models is also responsible for additional symmetries, as was already observed for the special case N = 2 [18]. Consider figure 4 with m = l (the same discussion holds for m = N − l). It is clear that the sum of the numbers appearing in the diamonds, on diagonals or anti-diagonals, is either 0, 2l or −2l, which is always even. This allows one to define conserved parity operators such as for example (−1) s∈red qs+ p∈pink qp , where the sum runs overs red stars and pink plaquettes forming a given diagonal. In order for this parity operator to be conserved when N is odd, one has to consider the charge and flux numbers q s and q p belonging to Z, introduced at the end of section 3.1 instead of q s and q p that belong to Z N . Of course, a dual discussion can be given by working on the dual lattice and in the eigenbasis of X l i Z m i operators. As in the Z 2 toric code in a transverse field, the conservation of these parity operators constrains the dynamics, and the model displays dimensional reduction. This dimensional reduction was originally discussed in the Xu-Moore model [35][36][37] which has the same spectrum as the Z 2 toric code in a transverse field. A detailed discussion about these issues can be found in [38].   Figure 6. Illustration (on the dual lattice defined in the right part of figure 5 ) of the action of the operators As and Bp on an eigenstate of X l i Z m i operators for all sites i. As in figure 4, a diamond on a star or a plaquette, with a value k = ±l, ±m denotes a multiplicative change by ω k of the corresponding eigenvalue. The behavior depends on the site's color (so on the star and plaquette types on the original lattice).

Model
In this section, we focus on the case N = 3 and we study the robustness of the Z 3 toric code with respect to simple perturbations. To this aim, and for the sake of simplicity, we consider the following Hamiltonian where we choose J = 1/3 in order to set the elementary excitation gap of the unperturbed Hamiltonian H TC to unity. In addition, we restrict our discussion to real parameters and set h 1,0 = h X , h 1,1 = h ⊥ , h 0,1 = h Z . For N = 2, this Hamiltonian (with the proper phase factor for H 1,1 discussed in section 3.3) corresponds to the Z 2 toric code in a uniform magnetic field studied in [19] so that this choice is well suited for a comparison between both systems. In the following, we adopt a language similar to that used for N = 2. Thus, H 1,0 and H 0,1 (H 1,1 ) will be considered as "parallel" ("transverse") perturbations.
The determination of the full three-dimensional phase diagram of Hamiltonian H(h X , h ⊥ , h Z ) is a difficult problem. Therefore, as was done initially for the Z 2 case [17,18], we will study parallel and transverse cases separately, but let us first discuss the methods used.

Methods
As already discussed, the system has to undergo phase transitions when perturbed with local operators. As in the perturbed Z 2 toric code, one expects first-and secondorder transitions in the phase diagram [19]. In order to analyze the breakdown of the topological phase, we combine the pCUT method and the variational iPEPS algorithm. This approach is motivated by the fact that the pCUT gives reliable estimates for second-order phase transitions whereas the iPEPS algorithm, as a variational tool, is especially sensitive to first-order transitions.

pCUT
The method of continuous unitary transformations has been introduced in reference [39] and general aspects of its perturbative variant pCUT can be found in reference [40]. In what follows, we focus on points that are specific to the application of the pCUT method to a topologically-ordered phase [17][18][19].
To apply the pCUT, it is essential that the unperturbed Hamiltonian (here H TC ) possesses an equidistant spectrum ‡ that is bounded from below [40]. These two constraints are satisfied in the Z 3 toric code as long as gaps of charges and fluxes are identical. In this case, one can interpret the toric code as a counting operator Q of charges and fluxes in the system ‡ We would like to stress that the spectrum of the Z N model studied in this paper is only equidistant for N = 2, 3, 4, so that one cannot use this machinery for other values of N . However, it could be applied to Kitaev's Z N toric code [3], whose spectrum is equidistant for all N .
where N s (N p ) denotes the total number of stars (plaquettes). Therefore the constant term represents the ground-state energy (remember that we set J = 1/3).
It is then possible to rewrite the local perturbations as n T n , where T n changes the particle number in the system by n, i.e., [Q, T n ] = n T n . The pCUT maps, order by order in the perturbation, the Hamiltonian H onto an effective Hamiltonian H eff , unitarily equivalent to H (same spectrum but different eigenstates), that reads as follows in the eigenbasis of the bare Hamiltonian H TC : As explained in [40], the coefficients C(m 1 , . . . , m k ) are model-independent rational numbers. An essential property of the effective Hamiltonian is that H eff , Q = 0. As a consequence, the number of quasiparticles (QPs) in the system, i.e., eigenstates of Q, is a good quantum number. In the perturbed toric code, QPs are dressed anyons adiabatically connected to the corresponding bare charges and fluxes.
To determine the zero-temperature phase diagram, we focus on the low-energy spectrum of H eff . Essentially, one must study H eff in the 0-QP subspace (to compute the ground-state energy) and in the 1-QP subspace (to obtain the low-energy gap ∆). We emphasize that computing the low-energy gap from the 1-QP subspace is meaningful as long as there are no bound states with lower energies. This is a working hypothesis that is crucial in what follows. As discussed in section 2.2.2, for open boundary conditions, the ground-state is non degenerate and the ground-state energy is nothing but the expectation value of H eff on the ground state of the bare Hamiltonian The structure of the 1-QP subspace is less trivial. Indeed, for N = 3 there are four different kinds of excitations : charges (q s = 1) and anti-charges (q s = −1 = 2 mod 3) living on stars, as well as fluxes (q p = 1) and anti-fluxes (q p = −1 = 2 mod 3) living on plaquettes. The associated subspaces are not connected by H eff because of the symmetry ensuring that the total charge and total flux are conserved (modulo N = 3 when working with the physical charge and flux). Thus H eff is a 1-QP hopping Hamiltonian in each of these sectors, and is therefore easily diagonalized, once the hopping amplitudes are determined. One then obtains dispersion relations for the four kinds of excitations, which only give two different energies because charges and anti-charges, as well as fluxes and anti-fluxes, play symmetric roles. Finally, one can compute the gap ∆ as the minimum over all momenta of these two energy bands.
Let us point out that the major challenge, in both the 0-QP and 1-QP sectors, lies in the computation of matrix elements of H eff . Indeed, one has to take into account the non-trivial braiding phases coming from virtual fluctuations of the excitations. Although the method yields results in the thermodynamical limit, the linked-cluster theorem (see for example [41,42] and references therein) allows one to compute the hopping amplitudes by considering finite-size clusters (albeit one needs a growing number of them when the order of the perturbation increases).
At the end of the day, one obtains a high-order series expansion of the groundstate energy E 0 and of the 1-QP gap ∆. The extrapolation of ∆ with standard resummation techniques (see e.g. [43]) allows a reliable determination of second-order phase transition points and thus of the boundaries of the topological phase (assuming that bound states of elementary QPs are not relevant and do not have an energy smaller than that of a single QP).
Unfortunately, as already stated, series expansions are not adapted to detect possible first-order phase transitions (in particular when having series in one phase only) so that one needs a complementary tool which we now describe.

iPEPS
The so-called iPEPS algorithm [44] is a variational method which, as such, is aimed at approximating the ground state of two-dimensional quantum lattice systems by employing a tensor-network approach. Details about this method have already been extensively discussed in the literature [44][45][46][47][48]. For completeness, though, we explain some of the basic features of the algorithm, focusing on those that are relevant for the study of the perturbed Z 3 toric code.
In the iPEPS algorithm, the quantum state |Ψ of the infinite square lattice is represented by a projected entangled pair state (PEPS) [44,49]. In the present problem, we choose a PEPS with four tensors denoted P , Q, R, and S per unit cell (see figure 7). Each of these tensors depends on O(d D 4 ) complex coefficients, where d = N = 3 is the dimension of the local Hilbert space at each site, and D is the so-called bond dimension of the PEPS. This bond dimension controls the maximum amount of entanglement carried by the PEPS wave function and, consequently, the accuracy of the ansatz. Following the discussion in [50] for N = 2, one can show that the ground state of the (non-perturbed) Z N toric code is a D = N -PEPS. Obviously, for J = 0 the fully-polarized ground state is also a (trivial) PEPS (D = 1). Consequently, at least in these two limiting cases, PEPS are exact ground states of the Hamiltonian.
In practice, for a given D, the goal is to find the coefficients of tensors P , Q, R and S that best approximate the ground state of H. These coefficients can be determined by an imaginary-time evolution driven by the Hamiltonian, since where |Ψ gs is the ground state of H and |Ψ 0 is any initial state that has a nonvanishing overlap with the ground state. The approximation of this evolution is performed in a similar way as explained, for instance, in [44,47].
(i) The whole evolution is splitted into small imaginary-time steps δτ by using a Suzuki-Trotter expansion of the evolution operator e −τ H . More precisely, writing the Hamiltonian as a sum of four-body terms h [ijkl] , we consider, at each time step, the action of the four-site operator (ii) At each imaginary-time step the state is approximated by some PEPS with the considered structure and bond dimension D. For instance, if at step τ we have a PEPS |Ψ(τ ) , then the evolved state | Ψ(τ + δτ ) ≡ g [ijkl] |Ψ(τ ) is also approximated by a new PEPS |Ψ(τ + δτ ) with the same structure. Practically, this approximation is achieved by minimizing the distance ||| Ψ(τ + δτ ) − |Ψ(τ + δτ ) || 2 with respect to the coefficients of tensors P , Q, R and S of the new PEPS.
In our case, we have carried out this minimization simultaneously over the four tensors by using a standard conjugate-gradient algorithm. Quite importantly, step (ii) as well as the evaluation of expectation values of local observables involves the contraction of an infinite two-dimensional tensor network. This contraction can be approximated by various schemes [44,47,51]. Here, we choose the Directional Corner Transfer Matrix Approach introduced in [47] that can be easily adapted to deal with different types of two-dimensional tensor networks [48], including those considered here. An important parameter in these manipulations is the so-called bond dimension of the environment χ which controls the accuracy of the approximations involved at this step.
Thus, according to the discussion above, there are four possible sources of error in the iPEPS algorithm. In this study, we have fixed D = 3, χ = 30, δτ = 10 −3 , and the aforementioned "four-site" unit cell. We checked that the increase of precision obtained by varying the values of the parameters is within the error bars obtained by the pCUT approach. Therefore, this choice of parameters turns out to be sufficient for our purposes.

pCUT+iPEPS
To determine the boundaries of the topological phase, we combine results from pCUT and iPEPS algorithm as explained below. To simplify the discussion, let us assume that we have only one control parameter h > 0.
Let us recall that a second-order transition is associated to the closure of the low-energy gap. Assuming that the relevant gap comes from the 1-QP sector, the critical point h c is then defined by ∆ (h c ) = 0. As already mentioned, this point can be efficiently computed with the pCUT method by extrapolating the high-order series expansion of ∆. Typically, with the maximum orders reached in this study, one can determine h c with a relative precision of the order 10 −2 -10 −3 . However, if some level crossing occurs, one faces a first-order transition that cannot be captured by the criterion ∆ = 0. This is why it is crucial to use the iPEPS algorithm to compute, variationally, the ground-state energy. Indeed, denoting by e pCUT 0 and e iPEPS 0 the ground-state energies calculated by both methods and assuming the existence of a point h where e iPEPS 0 < e pCUT 0 , two cases must be considered. If h > h c , it means that the iPEPS algorithm does not detect any level crossing for the ground state before the critical point, and hence a second order transition is likely taking place at h c . By contrast, if h < h c , then it means that a level crossing definitely occurs before the critical point h c and we conclude that a first-order transition takes place at h .
Of course, this reasoning would be exact if (i) one would have an infinite precision in both methods which is clearly not the case and (ii) no bound-state with lower energy exists (see discussion above). The accuracy of this combined pCUT+iPEPS approach is limited by several sources. First, the series expansion is performed up to a finite maximum order and the error of resummation schemes like the dlog-Padé extrapolation is hard to quantify. Second, the variational iPEPS algorithm is limited by the values of the parameters D, χ, δτ as well as the structure of the tensor network itself (see section 4.2.2). Additionally, it is numerically challenging to extract the global minimum of the variational ground-state energy. Reasonably, one can state that the combined pCUT+iPEPS approach works well as long as the error bars of both methods when determining the values h and h c are small compared to the difference |h − h c |. As explained above, in the present work, this relative error on h and h c is of the order 10 −2 -10 −3 .

Results
Let us now present our results concerning the perturbed Z 3 toric code H(h X , h ⊥ , h Z ) [see (22)] for several limiting cases. 4.3.1. The case h X = h ⊥ = 0 As already mentioned in section 3.2, the low-energy physics of the perturbed toric code H(0, 0, h Z ) corresponds to the N -state clock model in a transverse field (18). For N = 3, this model is also equivalent to the three-state Potts model in a transverse field. The extension of the topological phase can therefore be obtained by directly analyzing this model, which is simpler since one only has to consider charge degrees of freedom that live on stars of the square lattice. Indeed, remind that fluxes are frozen (conserved) and even absent in the low-energy sector.
As a first step, apart from the pCUT+iPEPS analysis, let us perform a standard mean-field calculation that, as we shall see, already captures the main qualitative aspects of the phase diagram although it relies on a trivial (non-entangled) variational state. To this aim, let us consider the following trial wave function where the coefficients a b,r , b b,r , and c b,r are chosen such that the local wave functions |ψ b,r s are normalized and minimize Ψ|H|Ψ . The introduction of different wave functions on red (r) and blue (b) stars is needed to accomodate with both ferromagnetic order and anti-ferromagnetic order (expected for h Z > 0 and h Z < 0 respectively). For the clock model, the sublattice magnetization is a proper order parameter. The topological phase of the perturbed Z 3 toric code corresponds to the disordered (symmetric) phase of the three-state clock model characterized by M b,r = 0. This phase is obtained for |h Z | J = 1/3. By contrast, an ordered (broken) phase with M b,r = 0 is expected for large perturbations |h Z | J = 1/3. Once again, let us emphasize that the sign of h Z is important for N = 3. Indeed, for large positive h Z , a polarized phase (with all stars in the same state) is stabilized, whereas for large negative h Z a "staggered" order arises.
Our results for the mean-field order parameter are summarized in the upper panel of figure 8. To go beyond this mean-field analysis and to obtain more quantitative results, let us now turn to the pCUT+iPEPS analysis. As the Hamiltonian (18) only contains one-site and two-site terms, the evaluation of effective matrix elements in the pCUT calculation can be performed using a full graph decomposition [43,52] that allowed us to reach order 11 (see Appendix A). The iPEPS algorithm also considerably benefits from the absence of four-site terms in H clock . As can be seen in Fig. 8, the PCUT+iPEPS method is in qualitative agreement with the mean-field treatment since it also predicts a first-order transition for h Z > 0 and a second-order transition for h Z < 0. However, for h Z > 0 where our results match those given in [53] for the Potts model (up to a rescaling), we find h Z 0.126 and h c Z 0.129. Thus, within the combined pCUT+iPEPS scheme, we are led to conclude that a (weakly) first-order phase transition occurs at h Z 0.126 in agreement with the results given in [28,46]. As can be seen in Fig. 8 (upper panel), the calculation of the order parameter using the iPEPS algorithm confirms a first-order behavior (jump of the magnetization). Note also that the relative difference between h Z and h 1,MF Z is about 12%. The situation is different for h Z < 0. In this case, we find h −0.204 and h c Z −0.195 (see table A1 for more details). Let us stress that we observed a very good agreement between iPEPS and PCUT calculations with a relative error between both ground-state energies smaller than 10 −4 for h c Z < h Z < 0. We therefore conclude that, within our scheme, a second-order phase transition occurs at h c Z −0.195. In this parameter region, the discrepancy between h c Z and h 2,MF Z = −1/8 is larger than 35%. To the best of our knowledge, this second-order phase transition has never been discussed in the literature and it is therefore very desirable to further characterize its universality class. Unfortunately, the iPEPS calculation of the critical exponent β associated to the order parameter is known to be specially sensitive to finite-D effects. As shown in the one-dimensional Ising model in a transverse field [54], this exponent also eventually reaches its mean-field value 1/2 for any finite D when approaching the critical point. Nevertheless, the exact value (1/8 in the latter problem) can be observed in a field range near the critical point whose size increases with D. In a two-dimensional system, it is very difficult to increase D in order to perform a similar study. For the problem at hand, our fixed bond dimension D = 3 only allowed us to see an exponent different from β MF = 1/2 in a very small region, and this was not sufficient to determine a reliable value.
Alternatively, dlog-Padé extrapolations of high-order series expansion of ∆ allows one to determine the exponent driving the closure of the gap at the critical point. More precisely, denoting z the dynamical exponent and ν the correlation length exponent, one has ∆ ∼ |h − h c | z ν for h near h c (see e.g. [55]). At order 11, we found z ν 0.71 but, as can be seen in table A1, it is clear that this value is poorly converged. However, to roughly estimate the error, one may draw a parallel with the N = 2 problem (i.e., the Ising model in a transverse field) for which series expansion of the gap have been computed up to order 13 in [56]. There, using the order 11 results, one gets an exponent z ν 0.645 which only differs by a few percent from the commonly accepted values z = 1 and ν = 0.630(1). Hopefully, for N = 3, we are also close from the asymptotic value but one clearly needs a more quantitative study to clarify the nature of this quantum phase transition. From that respect, Monte-Carlo simulations could provide valuable insights.

4.3.2.
The case h ⊥ = 0 Let us now turn to the case when both perturbations H 0,1 and H 1,0 are present so that neither charges nor fluxes are locally conserved anymore. Thus, one has to treat charges and fluxes at the same level and to carefully take their mutual statistics into account in the virtual braiding processes. This constraint strongly reduces the maximum order that can be reached in the pCUT calculation. Instead of order 11, we only computed e 0 and ∆ up to order 7 for this case (see Appendix A). The iPEPS calculation becomes more involved as well since one now has to deal with 4-body interactions (instead of 2-body interactions in the three-state clock model). In other words, our results are less accurate when a more complex perturbation is considered, as expected.
In figure 9 (left), we display the phase diagram obtained by combining pCUT and iPEPS algorithm following the same procedure as previously but for arbitrary directions in the (h X , h Z ) plane. As can be seen, the shape of the topological phase is symmetric under the exchange of h X ↔ h Z . This is due to the fact that the perturbation H 0,1 + H 1,0 respects the "charge-flux symmetry" present in H TC . The transition lines that mark the boundaries of this phase are found to be either first-order or second-order lines.
Let us start with a discussion of the (two symmetric) first-order (cyan) lines that are directly connected to the weakly first-order transition points (cyan diamonds) of the three-state clock model in a (positive) transverse field. Near these points, h c and h are found to be very close so that it is challenging to clearly decide whether the transition is first-or second-order. However, given the precision reached and the existence of some well-identified first-order points (not shown in figure 9) away from the cyan diamonds, it seems reasonable to argue that the two lines are likely first order.
The situation is different for the other part of the topological phase whose boundaries are connected to the second-order points (black diamonds) of the threestate clock model in a (negative) transverse field. Here, the pCUT+iPEPS approach is always clearly consistent with a second-order phase transition.
In this region, the intersection of the two second-order transition lines (green triangle) is reminiscent of the phase diagram of the perturbed Z 2 toric code in a parallel magnetic field where a multi-critical point was found [16,17]. In the latter case, this crossing point is also connected to a finite-length first-order line that lies outside the topological phase. For the present problem (N = 3), we have not performed a small-J perturbation theory likely to reveal a similar feature. Nevertheless, we computed the gap exponent z ν along the second-order (black) lines. As can be seen in figure 9 (right), the behavior of this exponent is rather flat (same values as for h X = 0 or h Z = 0) except at its extremities where it increases significantly. This may indicate a different universality class at the crossing points (green triangle and magenta squares) as also found in the N = 2 problem. Anyway, let us stress that there are some limitations of our approach for the current problem. First, the perturbative expansion at order 7 is not sufficient to determine the gap exponent z ν accurately. Second, the finite width in h c X where the gap exponent differs from the value at h c X = 0 is likely an artifact of the finite-order series. Indeed, we rather expect that all points except crossing points belong to the same universality class as the anti-ferromagnetic clock model in a transverse field but one cannot exclude other scenarios. Once again, to obtain more quantitative results, it would be very valuable to perform numerical  [4,2] dlogPade [6,4] 3d XY Figure 9. Left : zero-temperature phase diagram of H(h X , 0, h Z ) obtained by the pCUT+iPEPS approach. Cyan (black) lines denote first-order (secondorder) phase transitions. The first-order (second-order) phase transitions for the ferromagnetic (anti-ferromagnetic) three-state clock model in a transverse field are marked by cyan (black) diamonds. Magenta squares locate the crossings between first-and second-order transition lines while the crossing of the two second-order lines is taking place at the green triangle. Right : plot of the gap exponent z ν as a function of h c X along the horizontal second-order phase transition line (see left part). The black solid line is the result obtained from the order seven expression of the 1-QP gap. The black square is our best estimate of z ν for the simpler case of one parallel perturbation using the order eleven series expansion obtained for the three-state Potts model in a transverse field (see section 4.3.1). simulations of this model by means of alternative methods.

4.3.3.
The case h X = h Z = 0 To conclude this analysis, let us consider the case of a "transverse" perturbation that, as discussed in section 3.3 for general N , leads to a self-dual spectrum for H(0, h ⊥ , 0). Since J > 0 and N is odd, one must also have h ⊥ > 0. A convenient parametrization for this problem consists in setting J = 1 3 cos θ and h ⊥ = 1 3 sin θ with θ ∈ [0, π/2]. The unperturbed toric code corresponds to θ = 0 whereas for θ = π/2 the Hamiltonian is purely local. The self-dual point is located at θ = π/4 (J = h ⊥ ) so that the spectrum is symmetric under the transformation θ ↔ π/2 − θ. Thanks to self-duality, one can determine directly the nature of the transition by studying the singularity of the ground-state energy. Note that we could not use this criterion in previous cases since we did not have reliable informations outside the topological phase. As can be clearly seen in figure 10 (lower panel), the ground-state energy displays a kink at the self-dual point ( ∂e0 ∂θ | θ=π/4 is discontinuous) so that a first-order transition occurs there.
As explained in section 3.3, the parity conservation rules constrain the dynamics. In particular, they prevent a single QP to move in the presence of such a perturbation. Consequently the 1-QP energy level of the toric code (θ = 0) does not give rise to dispersive bands but is simply renormalized when h ⊥ = 0. The corresponding gap ∆ is shown in figure 10 (upper panel) and does not vanish at the self-dual point. However, this observation is compatible with the existence of a first-order transition that is due to level crossings which cannot be captured by analyzing low-energy levels. This situation is exactly the same as the one discussed in [18] for the Z 2 toric code in a transverse field. As in [18], 2-QP bound states are either pinned or mobile in one dimension only, while the simplest two-dimensional dispersing object is a 4-QP bound state.

Conclusion and perspectives
In this work, we have introduced exactly solvable models with Z N (Abelian) anyons that generalize Kitaev's famous toric code [3] or, equivalently, Wen's plaquette model [27]. Our main motivation was to probe the robustness of Z N >2 topologically-ordered phases following recent works on the case N = 2 [14][15][16][17][18][19]. This question is crucial since such phases are believed to be protected against external perturbations up to an order proportional to the typical system size. Besides, this robustness has been proven recently for any local perturbation [4,5]. However, as early noticed by Kitaev in his seminal paper [3] : "Of course, the perturbation should be small enough, or else a phase transition may occur." To investigate the limits of the standard aforementioned perturbative argument, we added local perturbations to this Z N toric code and we obtained several important results. First of all, for specific choices of the perturbation, the perturbed Z N Hamiltonian can be mapped onto the two-dimensional N -state quantum clock model in a transverse field. This mapping generalizes the correspondence between the Z 2 toric code in a parallel field and the transverse-field Ising model [14,15] to arbitrary N . For N = 3 and antiferromagnetic couplings we found evidence for a second-order phase transition but we have not been able to determine accurately its universality class. Let us simply note that our estimate of the critical exponent z ν is compatible with that of the three-dimensional classical XY model describing the three-dimensional threestate classical antiferromagnetic Potts model (see, e.g., [57]). Second, we have also shown that for N = 3, multi-critical points are likely present in the phase diagram of H(h X , 0, h Z ) [defined in (22)] although, here again, it would be important to analyze them with complementary tools. In particular, numerical simulations would be as valuable as for N = 2 [16]. In addition, as already observed for the N = 2 case [18], we have shown that self-dual Hamiltonians may arise provided the perturbating operators satisfy some constraints that we derived for arbitrary N . In this case, the self-duality is associated with a dimensional reduction and, for N = 3, we found that a first-order transition occurs at the self-dual point (as was also the case for N = 2).
From a methodological point of view, the combination of high-order perturbation theory (pCUT) and variational techniques (iPEPS) has been shown to be very efficient in a domain where few alternative approaches exist. From that respect, let us mention that improvement of the variational ansatz could certainly be achieved by taking into account the gauge symmetry in the tensor network as recently discussed in [20,58,59].
To conclude, we would like to mention related problems that would be worth investigating in order to deepen our understanding of topological phases' fragility. In the simplest case N = 2, it would already be worthwhile considering a perturbed toric code on a lattice which is not self-dual such as, for instance, the honeycomb lattice. Indeed, in this case, the role of charges and fluxes degrees of freedom cannot be exchanged by a simple lattice transformation and it is likely that such a model could lead to a non-trivial phase diagram in the presence of a uniform magnetic field. One may also directly perturb Wen's plaquette model with arbitrary local perturbations. Indeed, as shown in [60] for N = 2, although Wen's and Kitaev's models are (almost) the same in the absence of perturbation, they display very different properties in the presence of a uniform magnetic field. To go beyond N = 2, we believe that the study of the breakdown of Z N topological phases initiated here would greatly benefit from a large-N analysis that might be performed in a field-theoretical framework. Although it is clearly beyond the scope of this paper it may shed light on several issues as for instance the universality class that can be met at the boundaries of topological phases. Let us mention that recent studies suggested the possibility that conformal quantum critical points may describe these phase transitions (see for instance [61,62]). Finally, the most important step certainly consists in analyzing the robustness of non-Abelian topological phases [22,63] which are of direct interest for topological quantum computation but, undoubtedly, these quantum objects are much more tricky to handle. Since this model is self-dual, the corresponding quantities in the opposite limit h ⊥ J > 0 are obtained by changing θ into π 2 − θ.