Preparing topologically ordered states by Hamiltonian interpolation

We study the preparation of topologically ordered states by interpolating between an initial Hamiltonian with a unique product ground state and a Hamiltonian with a topologically degenerate ground state space. By simulating the dynamics for small systems, we numerically observe a certain stability of the prepared state as a function of the initial Hamiltonian. For small systems or long interpolation times, we argue that the resulting state can be identified by computing suitable effective Hamiltonians. For effective anyon models, this analysis singles out the relevant physical processes and extends the study of the splitting of the topological degeneracy by Bonderson. We illustrate our findings using Kitaev's Majorana chain, effective anyon chains, the toric code and Levin-Wen string-net models.


Introduction
Topologically ordered phases of matter have attracted significant interest in the field of quantum information, following the seminal work of Kitaev [Kit03]. From the viewpoint of quantum computing, one of their most attractive features is their ground space degeneracy: it provides a natural quantum error-correcting code for encoding and manipulating information. Remarkably, the ground space degeneracy is approximately preserved in the presence of weak static Hamiltonian perturbations [BHM10,BH11,MZ13]. This feature suppresses the uncontrolled accumulation of relative phases between code states, and thus helps to overcome decoherence. This is a necessary requirement for the realization of many-body quantum memories [DKLP02].
To use topologically ordered systems as quantum memories and for fault-tolerant quantum computation, concrete procedures for the preparation of specific ground states are required. Such mechanisms depend on the model Hamiltonian which is being realized as well as on the particular experimental realization. Early work [DKLP02] discussed the use of explicit unitary encoding circuits for the toric code. This consideration is natural for systems where we have full access to unitary gates over the underlying degrees of freedom. We may call this the bottom-up approach to quantum computing: here one proceeds by building and characterizing individual components before assembling them into larger structures. An example are arrays of superconducting qubits [BKM + 14, CGM + 14, CMS + 15]. Other proposed procedures for state preparation in this approach involve engineered dissipation [DKP14, BBK + 13] or measurementbased preparation [LMGH15]. However, achieving the control requirements for experimentally performing such procedures is quite challenging. They require either a) independently applying complex sequences of gates on each of the elementary constituents b) precisely engineering a dissipative evolution, or c) performing an extensive set of local measurements and associated non-local classical data processing to determine and execute a suitable unitary correction operation. Imperfections in the implementation of such protocols pose a severe problem, especially in cases where the preparation time is extensive [BHV06,KP14].
In fact, these procedures achieve more than is strictly necessary for quantum computation: any ground state can be prepared in this fashion. That is, they constitute encoders, realizing an isometry from a number of unencoded logical qubits to the ground space of the target Hamiltonian. We may ask if the task of preparing topologically ordered state becomes easier if the goal is to prepare specific states instead of encoding arbitrary states. In particular, we may ask this question in the top-down approach to quantum computing, where the quantum information is encoded in the ground space of a given condensed matter Hamiltonian. An example are Majorana wires [MZF + 12, NPDL + 14] or fractional quantum Hall substrates [VYPW11]. Indeed, a fairly standard approach to preparing ground states of a Hamiltonian is to cool the system by weakly coupling it with a thermal bath at a temperature significantly lower than the Hamiltonian gap. Under appropriate ergodicity conditions, this leads to convergence to a state mainly supported on the ground space. Unfortunately, when using natural equilibration processes, convergence may be slow, and the resulting prepared state is generally a (logical) mixed state unsuitable for computation.
A natural alternative method for preparing ground states of a given Hamiltonian is adiabatic evolution: here one initializes the system in an easy-to-prepare state (e.g., a product state), which is the unique ground state of a certain initial Hamiltonian (e.g., describing a uniform field). Subsequently, the Hamiltonian of the system is gradually changed (by tuning external control parameters in a time-dependent fashion) until the target Hamiltonian is reached. If this time-dependent change of the Hamiltonian is "slow enough", i.e., satisfies a certain adiabaticity condition (see Section 2), the state of the system will closely follow the trajectory of instantaneous ground states. The resulting state then is guaranteed to be mainly supported on the ground space of the target Hamiltonian, as desired.
Adiabatic preparation has some distinct advantages compared to e.g., encoding using a unitary circuit. For example, in contrast to the latter, adiabatic evolution guarantees that the final state is indeed a ground state of the actual Hamiltonian describing the system, independently of potential imperfections in the realization of the ideal Hamiltonians. In contrast, a unitary encoding circuit is designed to encode into the ground space of an ideal model Hamiltonian, and will therefore generally not prepare exact ground states of the actual physical system (which only approximate the model Hamiltonian). Such an encoding into the ideal ground space may lead to a negligible quantum memory time in the presence of an unknown perturbation [PKSC10]; this is because ideal and non-ideal (perturbed) ground states may differ significantly (this phenomenon is referred to as Anderson's orthogonality catastrophe [And67]). Adiabatic evolution, on the other hand, elegantly sidesteps these issues.
The fact that adiabatic evolution can follow the actual ground state of a system Hamiltonian makes it a natural candidate for achieving the task of topological code state preparation. An additional attractive feature is that its experimental requirements are rather modest: while some time-dependent control is required, this can be local, and additionally translation-invariant. Namely, the number of external control parameters required does not scale with the system size or code distance.

Summary and outlook
Motivated by these observations, we consider the general problem of preparing topologically ordered states by what we refer to as Hamiltonian interpolation. We will use this terminology instead of "adiabatic evolution" since in some cases, it makes sense to consider scenarios where adiabaticity guarantees cannot be given. For concreteness, we consider a time-dependent Hamiltonian H(t) which monotonically sweeps over the path i.e., we assume that the interpolation is linear in time and takes overall time 1 T . Guided by experimental considerations, we focus on the translation-invariant case: here the Hamiltonians H(t) are translation-invariant throughout the evolution. More precisely, we consider the process of interpolating between a Hamiltonian H triv with unique ground state Ψ(0) = ϕ ⊗L and a Hamiltonian H top with topologically degenerate ground space (which is separated from the remainder of the spectrum by a constant gap): the state Ψ(t) of the system at time t ∈ [0, T ] satisfies the equation of motion Generally, we consider families of Hamiltonians (or models) parametrized by a system size L; throughout, we will assume that L is the number of single particles, e.g., the number of qubits (or sites) in a lattice with Hilbert space H = (C 2 ) ⊗L . The dimension of the ground space of H top will be assumed to be independent of the system size.
Our goal is to characterize the set of states which are preparable by such Hamiltonian interpolations starting from various product states, i.e., by choosing different initial Hamiltonians H triv . To each choice Ψ(0) = ϕ ⊗L of product state we associate a normalized initial trivial Hamiltonian H triv := − j P (j) ϕ which fully specifies the interpolating path of Eq. (1), with P (j) ϕ = |ϕ ϕ| being the single particle projector onto the state ϕ at site j.
In the limit T → ∞, one may think of this procedure as associating an encoded (logical) state ι(ϕ) to any single-particle state ϕ. However, some caveats are in order: first, the global phase of the state ι(ϕ) cannot be defined in a consistent manner in the limit T → ∞, and is therefore not fixed. Second, the final state in the evolution (2) does not need to be supported entirely on the ground space of H top because of non-adiaticity errors, i.e., it is not a logical (encoded) state itself. To obtain a logical state, we should think of ι(ϕ) as the final state projected onto the ground space of H top . Up to these caveats, our goal is essentially to characterize the image of the association ι : ϕ → ι(ϕ), as well as its continuity properties. We will also define an analogous map ι T associated to fixed evolution time T and study it numerically by simulating the corresponding Schrödinger equation (2) on a classical computer.
While there is a priori no obvious relationship between the final states ι T (ϕ), ι T (ϕ ) resulting from different initial (product) states ϕ ⊗L , ϕ ⊗L , we numerically find that the image of ι T is concentrated around a particular discrete family of encoded states. In particular, we observe for small system sizes that the preparation enjoys a certain stability property: variations in the initial Hamiltonian do not significantly affect the final state. We support this through analytic arguments, computing effective Hamiltonians associated to perturbations around H top which address the large T limit. This also allows us to provide a partial prediction of which states ι(ϕ) may be obtained through such a preparation process. We find that under certain general conditions, ι(ϕ) belongs to a certain finite family of preferred states which depend on the final Hamiltonian H top . As we will argue, there is a natural relation between the corresponding states ι(ϕ) for different system sizes: they encode the same logical state if corresponding logical operators are chosen (amounting to a choice of basis of the ground space).
Characterizing the set {ι(ϕ)} ϕ of states preparable using this kind of Hamiltonian interpolation is important for quantum computation because certain encoded states (referred to as "magic states") can be used as a resource for universal computation [BK05]. Our work provides insight into this question for 'small' systems, which we deem experimentally relevant. Indeed, there is a promising degree of robustness for the Hamiltonian interpolation to prepare certain (stabilizer) states. However, a similar preparation of magic states seems to require imposing additional symmetries which will in general not be robust. We exemplify our considerations using various concrete models, including Kitaev's Majorana chain [Kit01] (for which we can provide an exact solution), effective anyon chains (related to the so-called golden chain [FTL + 07] and the description used by Bonderson [Bon09]), as well as the toric code [Kit03] and Levin-Wen string-net models [LW05] (for which we simulate the time-evolution for small systems, for both the doubled semion and the doubled Fibonacci model).

Prior work
The problem of preparing topologically ordered states by adiabatic interpolation has been considered prior to our work by Hamma and Lidar [HL08]. Indeed, their contribution is one of the main motivations for our study. They study an adiabatic evolution where a Hamiltonian having a trivial product ground state is interpolated into a toric code Hamiltonian having a four-fold degenerate ground state space. They found that while the gap for such an evolution must forcibly close, this may happen through second order phase transitions. Correspondingly, the closing of the gap is only polynomial in the system size. This allows an efficient polynomialtime Hamiltonian interpolation to succeed at accurately preparing certain ground states. We revisit this case in Section 2.1 and give further examples of this phenomenon. The authors of [HZHL08] also observed the stability of the encoded states with respect to perturbations in the preparation process.
Bonderson [Bon09] considered the problem of characterizing the lowest order degeneracy splitting in topologically ordered models. Degeneracy lifting can be associated to tunneling of anyonic charges, part of which may be predicted by the universal algebraic structure of the anyon model. Our conclusions associated to Sections 5 and 6 can be seen as supporting this perspective.

Beyond small systems
In general, the case of larger systems (i.e., the thermodynamic limit) requires a detailed understanding of the quantum phase transitions [Sac11] occurring when interpolating between H triv and H top . Taking the thermodynamic limit while making T scale as a polynomial of the system size raises a number of subtle points. A major technical difficulty is that existing adiabatic theorems do not apply, since at the phase transition gaps associated to either of the relevant phases close. This is alleviated by scaling the interpolation time T with the system size and splitting the adiabatic evolution into two regimes, the second of which can be treated using degenerate adiabatic perturbation theory [RO10,RO12,RO14]. However, such a methodology still does not yield complete information about the dynamical effects of crossing a phase boundary.
More generally, it is natural to conjecture that interpolation between different phases yields only a discrete number of distinct states corresponding to a discrete set of continuous phase transitions in the thermodynamic limit. Such a conjecture links the problem of Hamiltonian interpolation to that of classifying phase transitions between topological phases. It can be motivated by the fact that only a discrete set of possible condensate-induced continuous phase transitions is predicted to exist in the thermodynamic limit [BS09,BSS11].

Adiabaticity and ground states
The first basic question arising in this context is whether the evolution (2) yields a state Ψ(T ) close to the ground space of H top . The adiabatic theorem in its multiple forms (see e.g., [Teu03]) provides sufficient conditions for this to hold: These theorems guarantee that given a Hamiltonian path {H(t)} 0≤t≤T satisfying certain smoothness and gap assumptions, initial eigenstates evolve into approximate instantaneous eigenstates under an evolution of the form (2). The latter assumptions are usually of the following kind: (i) Uniform gap: There is a uniform lower bound ∆(t) ≥ ∆ > 0 on the spectral gap of H(t) for all t ∈ [0, T ]. The relevant spectral gap ∆(t) is the energy difference between the ground space P 0 (t)H of the instantaneous Hamiltonian H(t) and the rest of its spectrum.
Here and below, we denote by P 0 (t) the spectral projection onto the ground space 2 of H(t).
(ii) Smoothness: There are constants c 1 , . . . , c M such that the M first derivatives of H(t) are uniformly bounded in operator norm, i.e., for all j = 1, . . . , M , we have The simplest version of such a theorem is: Theorem 2.1. Given a state Ψ(0) such that P 0 (0)Ψ(0) = Ψ(0) and a uniformly gapped Hamiltonian path H(t) for t ∈ [0, T ] given by Eq.
(1), the state Ψ(T ) resulting from the evolution (2) satisfies In other words, in the adiabatic limit of large times T , the state Ψ(T ) belongs to the instantaneous eigenspace P 0 (T )H and its distance from the eigenspace is O(1/T ).
This version is sufficient to support our analytical conclusions qualitatively. For a quantitative analysis of non-adiabaticity errors, we perform numerical simulations. Improved versions of the adiabatic theorem (see [GMC15,LRH09]) provide tighter analytical error estimates for general interpolation schedules at the cost of involving higher order derivatives of the Hamiltonian path H(t) (see Eq.'(3)), but do not change our main conclusions.
Several facts prevent us from directly applying such an adiabatic theorem to our evolution (1) under consideration.
Topological ground space degeneracy. Most notably, the gap assumption (i) is not satisfied if we study ground spaces: we generally consider the case where H(0) = H triv has a unique ground state, whereas the final Hamiltonian H(T ) = H top is topologically ordered and has a degenerate ground space (in fact, this degeneracy is exact and independent of the system size for the models we consider). This means that if P 0 (t) is the projection onto the ground space of H(t), there is no uniform lower bound on the gap ∆(t).
We will address this issue by restricting our attention to times t ∈ [0, κT ], where κ ≈ 1 is chosen such that H(κT ) has a non-vanishing gap but still is "inside the topological phase". We will illustrate in specific examples how Ψ(T ) can indeed be recovered by taking the limit κ → 1.
We emphasize that the expression "inside the phase" is physically not well-defined at this point since we are considering a Hamiltonian of a fixed size. Computationally, we take it to mean that the Hamiltonian can be analyzed by a convergent perturbation theory expansion starting from the unperturbed Hamiltonian H top . The resulting lifting of the ground space degeneracy of H top will be discussed in more detail in Section 3. Dependence on the system size. A second potential obstacle for the use of the adiabatic theorem is the dependence on the system size L (where e.g., L is the number of qubits). This dependence enters in the operator norms (3), which are extensive in L -this would lead to polynomial dependence of T on L even if e.g., the gap were constant (uniformly bounded).
More importantly, the system size enters in the gap ∆(t): in the topological phase, the gap (i.e., the splitting of the topological degeneracy of H top ) is exponentially small in L for constant-strength local perturbations to H top , as shown for the models considered here by Bravyi, Hastings and Michalakis [BHM10]. Thus a naïve application of the adiabatic theorem only yields a guarantee on the ground space overlap of the final state if the evolution time is exponentially large in L. This is clearly undesirable for large systems; one may try to prepare systems faster (i.e., more efficiently) but would need alternate arguments to ensure that the final state indeed belongs to the ground space of H top .
For these reasons, we restrict our attention to the following two special cases of the Hamiltonian interpolation (1): • Symmetry-protected preparation: if there is a set of observables commuting with both H triv and H top , these will represent conserved quantities throughout the Hamiltonian interpolation. If the initial state is an eigenstate of such observables, one may restrict the Hilbert space to the relevant eigenvalue, possibly resolving the topological degeneracy and guaranteeing a uniform gap. This observation was first used in [HL08] in the context of the toric code: for this model, such a restriction allows mapping the problem to a transverse field Ising model, where the gap closes polynomialy with the system size. We identify important cases satisfying this condition. While this provides the most robust preparation scheme, the resulting encoded states are somewhat restricted (see Section 2.1).
• Small systems: For systems of relatively small (constant) size L , the adiabatic theorem can be applied as all involved quantities are essentially constant. In other words, although 'long' interpolation times are needed to reach ground states of H top (indeed, these may depend exponentially on L), these may still be reasonable experimentally. The consideration of small system is motivated by current experimental efforts to realize surface codes [KBF + 15]: they are usually restricted to a small number of qubits, and this is the scenario we are considering here (see Section 2.2).
Obtaining a detailed understanding of the general large L limiting behaviour (i.e., the thermodynamic limit) of the interpolation process (1) is beyond the scope of this work.

Symmetry-protected preparation
Under particular circumstances, the existence of conserved quantities permits applying the adiabatic theorem while evading the technical obstacle posed by a vanishing gap in the context of topological order. Such a case was considered by Hamma and Lidar [HL08], who showed that 8 certain ground states of the toric code can be prepared efficiently. We can formalize sufficient conditions in the following general way (which then is applicable to a variety of models, as we discuss below).
(iii) The final ground space has support on QP 0 (T )H = 0 Then QΨ(t) = Ψ(t), and the adiabatic theorem can be applied with lower bound ∆ on the gap, The proof of this statement is a straightforward application of the adiabatic theorem (Theorem 2.1) to the Hamiltonians QH triv and QH top in the restricted subspace QH. In the following sections, we will apply Observation 2.2 to various systems. It not only guarantees that the ground space is reached, but also gives us information about the specific state prepared in a degenerate ground space.
As an example of the situation discussed in Observation 2.2, we discuss the case of fermionic parity conservation in Section 4. This symmetry is naturally present in fermionic systems. We expect our arguments to extend to more general topologically ordered Hamiltonians with additional symmetries. It is well-known that imposing global symmetries on top of topological Hamiltonians provides interesting classes of systems. Such symmetries can exchange anyonic excitations, and their classification as well as the construction of associated defect lines in topological Hamiltonians is a topic of ongoing research [BSW11,KK12,BJQ13]. The latter problem is intimately related to the realization (see e.g., [BMD09,Bom15]) of transversal logical gates, which leads to similar classification problems [BK13, BBK + 14, Yos15b,Yos15a]. Thus we expect that there is a close connection between adiabatically preparable states and transversally implementable logical gates. Indeed, a starting point for establishing such a connection could be the consideration of interpolation processes respecting symmetries realized by transversal logical gates.
For later reference, we also briefly discuss a situation involving conserved quantities which -in contrast to Observation 2.2 -project onto excited states of the final Hamiltonian. In this case, starting with certain eigenstates of the corresponding symmetry operator Q, the ground space cannot be reached: Observation 2.3. Assume that Q, H triv , H top , Ψ(0) obey properties (i),(ii) and (iv) of Observation 2.2. If the ground space P 0 (T )H of H top satisfies QP 0 (T )H = 0 (i.e., is orthogonal to the image of Q), then the Hamiltonian interpolation cannot reach the ground space of H top , i.e., Ψ(T ), P 0 (T )Ψ(T ) = Ω(1).

9
The proof of this observation is trivial since Q is a conserved quantity of the Schrödinger evolution. Physically, the assumptions imply the occurrence of a level-crossing where the energy gap exactly vanishes and eigenvalue of Q restricted to the ground space changes. We will encounter this scenario in the case of the toric code on a honeycomb lattice, see Section 7.3.

Small-system case
In a more general scenario, there may not be a conserved quantity as in Observation 2.2. Even assuming that the ground space is reached by the interpolation process (1), it is a priori unclear which of the ground states is prepared. Here we address this question.
As remarked earlier, we focus on systems of a constant size L, and assume that the preparation time T is large compared to L. Generically, the Hamiltonians H(t) are then non-degenerate (except at the endpoint, t ≈ T , where H(t) approaches H top ). Without fine tuning, we may expect that there are no exact level crossings in the spectrum of H(t) along the path t → H(t) (say for some times t ∈ [0, κT ], κ ≈ 1). For sufficiently large overall evolution times T , we may apply the adiabatic theorem to conclude that the state of the system follows the (unique) instantaneous ground state (up to a constant error). Since our focus is on small systems, we will henceforth assume that this is indeed the case, and summarily refer to this as the adiabaticity assumption. Again, we emphasize that this is a priori only reasonable for small systems.
Under the adiabaticity assumption, we can conclude that the prepared state Ψ(T ) roughly coincides with the state obtained by computing the (unique) ground state ψ κ of H(κT ), and taking the limit κ → 1. In what follows, we adopt this computational prescription for identifying prepared states. Indeed, this approach yields states that match our numerical simulation, and provides the correct answer for certain exactly solvable cases. Furthermore, the computation of the states ψ κ (in the limit κ → 1) also clarifies the physical mechanisms responsible for the observed stability property of preparation: we can relate the prepared states to certain linear combination of string-operators (Wilson-loops), whose coefficients depend on the geometry (length) of these loops, as well as the amplitudes of certain local particle creation/annihilation and tunneling processes.
Since H(κT ) for κ ≈ 1 is close to the topologically ordered Hamiltonian H top , it is natural to use ground states (or logical operators) of the latter as a reference to express the instantaneous states ψ κ . Indeed, the problem essentially reduces to a system described by H top , with an additional perturbation given by a scalar multiple of H triv . Such a local perturbation generically splits the topological degeneracy of the ground space. The basic mechanism responsible for this splitting for topologically ordered systems has been investigated by Bonderson [Bon09], who quantified the degeneracy splitting in terms of local anyon-processes. We seek to identify low-energy ground states: this amounts to considering the effective low-energy dynamics (see Section 3). This will provide valuable information concerning the set {ι(ϕ)}.

Effective Hamiltonians
As discussed in Section 2.2, for small systems (and sufficiently large times T ), the state Ψ(κT ) in the interpolation process (1) should coincide with the ground state of the instantaneous Hamiltonian H(κT ). For κ ≈ 1, the latter is a perturbed version of the Hamiltonian H top , where the perturbation is a scalar multiple of H triv . That is, up to rescaling by an overall constant, we are concerned with a Hamiltonian of the form where H 0 = H top is the target Hamitonian and V = H triv is the perturbation. To compute the ground state of a Hamiltonian of the form (4), we use effective Hamiltonians. These provide a description of the system in terms of effective low-energy degrees of freedom.

Low-energy degrees of freedom
Let us denote by P 0 the projection onto the degenerate ground space of H 0 . Since H 0 is assumed to have a constant gap, a perturbation of the form (4) effectively preserves the lowenergy subspace P 0 H for small > 0, and generates a dynamics on this subspace according to an effective Hamiltonian H eff ( ). We will discuss natural definitions of this effective Hamiltonian in Section 3.3. For the purpose of this section, it suffices to mention that it is entirely supported on the ground space of H 0 , i.e., H eff ( ) = P 0 H eff ( )P 0 . As such, it has spectral decomposition where E eff 0 < E eff 1 < . . ., and where Π eff k ( ) = Π eff k ( )P 0 are commuting projections onto subspaces of the ground space P 0 H of H 0 . (Generally, we expect H eff ( ) to be non-degenerate such that K = dim P 0 H.) In particular, the effective Hamiltonian (5) gives rise to an orthogonal decomposition of the ground space P 0 H by projections {Π eff k ( )} K−1 k=0 . States in Π eff 0 ( )H are distinguished by having minimal energy. We can take the limiting projections as the perturbation strength goes to 0, setting In particular, the effective Hamiltonian H eff ( ) has ground space Π eff 0 (0)H in the limit → 0. Studying H eff ( ), and, in particular, the space Π eff 0 (0)H appears to be of independent interest, as it determines how perturbations affect the topologically ordered ground space beyond spectral considerations as in [Bon09].

Hamiltonian interpolation and effective Hamiltonians
The connection to the interpolation process (1) is then given by the following conjecture. It is motivated by the discussion in Section 2.2 and deals with the case where there are no conserved quantities (unlike, e.g., in the case of the Majorana chain, as discussed in Section 4).
Conjecture 1. Under suitable adiabaticity assumptions (see Section 2.2) the projection of the final state Ψ(T ) onto the ground space of H top belongs to Π eff 0 (0)H (up to negligible errors 3 ), i.e., it is a ground state of the effective Hamiltonian H eff ( ) in the limit → 0.
In addition to the arguments in Section 2.2, we provide evidence for this conjecture by explicit examples, where we illustrate how Π eff 0 (0)H can be computed analytically. We also verify that Conjecture 1 correctly determines the final states by numerically studying the evolution (1).
We remark that the statement of Conjecture 1 severly constrains the states that can be prepared by Hamiltonian interpolation in the large T limit: we will argue that the space Π eff 0 (0)H has a certain robustness with respect to the choice of the initial Hamiltonian H triv . In fact, the space Π eff 0 (0)H is typically 1-dimensional and spanned by a single vector ϕ 0 . Furthermore, this vector ϕ 0 typically belongs to a finite family A ⊂ P 0 H of states defined solely by H top . In particular, under Conjecture 1, the dependence of the final state Ψ(T ) on the Hamiltonian H triv is very limited: the choice of H triv only determines which of the states in A is prepared. We numerically verify that the resulting target states Ψ(T ) indeed belong to the finite family A of states obtained analytically.

Perturbative effective Hamiltonians
As discussed in Section 3.2, we obtain distinguished final ground states by computation of suitable effective Hamiltonians H eff ( ), approximating the action of H 0 + V on the ground space P 0 H of H 0 . In many cases of interest, computing this effective Hamiltonian (whose definition for the Schrieffer-Wolff-case we present in Appendix A.1) exactly is infeasible (The effective Hamiltonian for the Majorana chain (see Section 4) is an exception.).
Instead, we seek a perturbative expansion in terms of powers of the perturbation strength . This is particularly natural as we are interested in the limit → 0 anyway (see Conjecture 1). Furthermore, it turns out that such perturbative expansions provide insight into the physical mechanisms underlying the 'selection' of particular ground states. We remark that there are several different methods for obtaining low-energy effective Hamiltonians. The Schrieffer-Wolff method [SW66, BDL11] provides a unitary U such that H eff = U (H 0 + V )U † preserves P 0 H and can be regarded as an effective Hamiltonian. One systematically obtains a series expansion for the anti-Hermitian generator S of U = e S ; this then naturally gives rise to an order-by-order expansion of the effective Hamiltonian, where P 0 is the projection onto the ground space P 0 H of H 0 (explicit expressions are given in Appendix A.2).
Using the Schrieffer-Wolff method has several distinct advantages, including the fact that (i) the resulting effective Hamiltonian H eff , as well as the terms H (n) eff are Hermitian, and hence have a clear physical interpretation. This is not the case e.g., for the Bloch expansion [Blo58].
(ii) There is no need to address certain self-consistency conditions arising e.g., when using the Dyson equation and corresponding self-energy methods [ABD75,FW03] We point out that the series resulting by taking the limit n → ∞ in (6) has the usual convergence issues encountered in many-body physics: convergence is guaranteed only if V ≤ ∆, where ∆ is the gap of H 0 . For a many-body system with extensive Hilbert space (e.g., L spins), the norm V = Ω(L) is extensive while the gap ∆ = O(1) is constant, leading to convergence only in a regime where = O(1/L). In this respect, the Schrieffer-Wolff method does not provide direct advantages compared to other methods. As we are considering the limit → 0, this is not an issue (also, for small systems as those considered in our numerics, we do not have such issues either). We point out, however, that the results obtained by Bravyi et al. [BDL11] suggest that considering partial sums of the form (6) is meaningful even in cases in which the usual convergence guarantees are not given: indeed, [BDL11,Theorem 3] shows that the ground state energies of H eff are essentially local (for low orders n) when the method is applied to certain many-body systems, see [BDL11]. We will not need the corresponding results here, however.
Unfortunately, computing the Schrieffer-Wolff Hamiltonian H (n) eff generally involves a large amount of combinatorics (see [BDL11] for a diagrammatic formalism for this purpose). In this respect, other methods may appear to be somewhat more accessible. Let us mention in particular the method involving the Dyson equation (and the so-called 'self-energy' operator), which was used e.g., in [Kit06, Section 5.1] to compute 4-th order effective Hamiltonians. This leads to remarkably simple expressions of the form for the n-th order term effective Hamiltonian, where G = G(E 0 ) is the resolvent operator evaluated at the ground state energy E 0 of H 0 . In general, though, the expression (7) only coincides with the Schrieffer-Wolff-method (that is, (6)) up to the lowest non-trivial order.

Perturbative effective Hamiltonians for topological order
Here we identify simple conditions under which the Schrieffer-Wolff Hamiltonian of lowest nontrivial order has the simple structure (7). We will see that these conditions are satisfied for the systems we are interested in. In other words, for our purposes, the self-energy methods and the Schrieffer-Wolff method are equivalent. While establishing this statement (see Theorem 3.2 below) requires some work, this result vastly simplifies the subsequent analysis of concrete systems. The condition we need is closely related to quantum error correction [KL97]. In fact, this condition has been identified as one of the requirements for topological quantum order (TQO-1) in Ref. [BHM10]. To motivate it, consider the case where P 0 H is an error-correcting code of distance L. Then all operators T acting on less than L particles 4 have trivial action on the code space, i.e., for such T , the operator P 0 T P 0 is proportional to P 0 (which we will write as P 0 T P 0 ∈ CP 0 ). In particular, this means that if V is a Hermitian linear combination of singleparticle operators, then P 0 V n P 0 ∈ CP 0 for all n < L. The condition we need is a refinement of this error-correction criterion that incorporates energies (using the resolvent): Definition 3.1. We say that the pair (H 0 , V ) satisfies the topological order condition with parameter L if L is the smallest interger such that for all n < L, we have Here P 0 is the ground space projection of H 0 , Q 0 = I −P 0 is the projection onto the orthogonal complement, and G = G(E 0 ) is the resolvent (8) (supported on Q 0 H).
We remark that this definition is easily verified in the systems we consider: if excitations in the system are local, the resolvent operators and projection in a product of the form (9) can be replaced by local operators, and condition (9) essentially reduces to a standard error correction condition for operators with local support.
Assuming this definition, we then have the following result: Theorem 3.2. Suppose that (H 0 , V ) satisfies the topological order condition with parameter L.
Then the n-th order Schrieffer-Wolff effective Hamiltonian satisfies i.e., the effective Hamiltonian is trivial for these orders, and We give the proof of this statement in Appendix A.

The Majorana chain
In this section, we apply our general results to Kitaev's Majorana chain. We describe the model in Section 4.1. In Section 4.2, we argue that the interpolation process (2) is an instance of symmetry-protected preparation; this allows us to identify the resulting final state. We also observe that the effective Hamiltonian is essentially given by a 'string'-operator F , which happens to be the fermionic parity operator in this case. That is, up to a global energy shift, we have for a certain constant f depending on the choice of perturbation.
4 By particle we mean a physical constituent qubit or qudit degree of freedom.
Without loss of generality, we have chosen the normalization such that elementary excitations have unit energy. The Hamiltonian has a two-fold degenerate ground space. The Majorana operators c 1 and c 2L correspond to a complex boundary mode, and combine to form a Dirac fermion which commutes with the Hamiltonian. The operator a † a hence provides a natural occupation number basis {|g σ } σ∈{0,1} for the ground space P 0 H defined (up to arbitrary phases) by As a side remark, note that the states |g 0 and |g 1 cannot be used directly to encode a qubit. This is because they have even and odd fermionic parity, respectively, and thus belong to different superselection sectors. In other words, coherent superposition between different parity sectors are nonphysical. This issue can be circumvented by using another fermion or a second chain, see [BK12]. Since the conclusions of the following discussion will be unchanged, we will neglect this detail for simplicity. We remark that the Hamiltonian H top of Eq. (10) belongs to a one-parameter family of extensively studied and well-understood quantum spin Hamiltonians. Indeed, the Jordan-Wigner transform of the Hamiltonian (with g ∈ R an arbitrary parameter) is the transverse field Ising model where X j and Z j are the spin 1/2 Pauli matrices acting on qubit j, j = 1, . . . , L. This transformation allows analytically calculating the complete spectrum of the translation invariant chain for both periodic and open boundary conditions [Pfe70]. The Hamiltonian H I,g has a quantum phase transition at g = 1, for which the lowest energy modes in the periodic chain have an energy scaling as 1/L. The open boundary case has been popularized by Kitaev as the Majorana chain and has a unique low energy mode a (see Eq. (11)) which has zero energy for g = 0 and for finite 0 < g < 1, becomes a dressed mode with exponentially small energy (in L) and which is exponentially localized at the boundaries.

State preparation by interpolation
The second term in (12) may be taken to be the initial Hamiltonian H triv for the interpolation process. More generally, to prepare ground states of H top , we may assume that our initial Hamiltonian is a quadratic Hamiltonian with a unique ground state. That is, H triv is of the form where V is a real antisymmetric 2L × 2L matrix. We will assume that it is bounded and local (with range r) in the sense that where · denotes the operator norm. As shown in [BK12, Theorem 1], the Hamiltonian H top + H triv has two lowest energy states with exponentially small energy difference, and this lowestenergy space remains separated from the rest of the spectrum by a constant gap for a fixed (constant) perturbation strength > 0. Estimates on the gap along the complete path H(t) are, to the best of our knowledge, not known in this more general situation. Let us assume that Ψ(0) is the unique ground state of H triv and consider the linear interpolation (2). The corresponding process is an instance of the symmetry-protected preparation, i.e., Observation 2.2 applies in this case. Indeed, the fermionic parity operator commutes with both H triv and H top . Therefore, the initial ground state Ψ(0) lies either in the even-parity sector, i.e., F Ψ(0) = Ψ(0), or in the odd-parity sector (F Ψ(0) = −Ψ(0)). (Even parity is usually assumed by convention, since the fermionic normal modes used to describe the system are chosen to have positive energy.) In any case, the ±1 eigenvalue of the initial ground state with respect to F will persist throughout the full interpolation. This fixes the final state: Lemma 4.1. Under suitable adiabaticity assumptions (see Observation 2.2), the resulting state in the evolution (2) is (up to a phase) given by the ground state |g 0 or |g 1 , depending on whether the initial ground state Ψ(0) lies in the even-or odd-parity sector.
In particular, if H triv = − gi 2 L j=1 c 2j−1 c 2j is given by the second term in (12), we can apply the results of [Pfe70]: the gap at the phase transition is associated with the lowest energy mode (which is not protected by symmetry) and is given by λ 2 (H I,g=1 ) = 2 sin [π/(2L + 1)]. In other words, it is linearly decreasing in the system size L. Therefore, the total evolution time T only needs to grow polynomially in the system size L for Hamiltonian interpolation to accurately follow the ground state space at the phase transition. We conclude that translation-invariant Hamiltonian interpolation allows preparing the state |g 0 in a time T polynomial in the system size L and the desired approximation accuracy.
To achieve efficient preparation through Hamiltonian interpolation, one issue that must be taken into account is the effect of disorder (possibly in the form of a random site-dependent chemical potential). In the case where the system is already in the topologically ordered phase, a small amount of Hamiltonian disorder can enhance the zero temperature memory time of the Majorana chain Hamiltonian [BK12]. This 1D Anderson localization effect [And58], while boosting memory times, was also found to hinder the convergence to the topological ground space through Hamiltonian interpolation. Indeed, in [CFS07] it was found that the residual energy density [E res (T )/L] av ∝ 1/ ln 3.4 (T ) averaged over disorder realizations decreases only polylogarithmically with the Hamiltonian interpolation time. Such a slow convergence of the energy density indicates that in the presence of disorder, the time T required to accurately reach the ground space scales exponentially with the system size L. For this reason, translationinvariance (i.e., no disorder) is required for an efficient preparation, and this may be challenging in practice.
We emphasize that according to Lemma 4.1, the prepared state is largely independent of the choice of the initial Hamiltonian H triv (amounting to a different choice of V): we do not obtain a continuum of final states. As we will see below, this stability property appears in a similar form in other models. The parity operator (13), which should be thought of as a string-operator connecting the two ends of the wire, plays a particular role -it is essentially the effective Hamiltonian which determines the prepared ground state.
Indeed, the Schrieffer-Wolff-effective Hamiltonian can be computed exactly in this case, yielding where E 0 ( ) is the ground state energy of H top + H triv , and ∆( ) Expression (14) can be computed based on the variational expression (55) for the Schrieffer-Wolff transformation, using the fact that the ground space is two-dimensional and spanned by two states belonging to the even-and odd-parity sector, respectively. Note that the form (14) can also be deduced (without the exact constants) from the easily verified fact (see e.g., Eq. (54)) that the Schrieffer-Wolff unitary U commutes with the fermionic parity operator F , and thus the same is true for H eff ( ). This expression illustrates that Conjecture 1 does not directly apply in the context of preserved quantities, as explained in Section 3.2: rather, it is necessary to know the parity of the initial state Ψ(0) to identify the resulting final state Ψ(T ) in the interpolation process.

General anyon chains
In this section, we generalize the considerations related to the Majorana chain to more general anyonic systems. Specifically, we consider a 1-dimensional lattice of anyons with periodic boundary conditions. This choice retains many features from the Majorana chain such as locally conserved charges and topological degeneracy yet further elucidates some of the general properties involved in the perturbative lifting of the topological degeneracy. In Section 5.1, we review the description of effective models for topologically ordered systems. A key feature of these models is the existence of a family {F a } a of string-operators indexed by particle labels. Physically, the operators F a correspond to the process of creating a particle-antiparticle pair (a,ā), tunneling along the 1-dimensional (periodic) lattice, and subsequent fusion of the pair to the vacuum (see Section 5.1.6). These operators play a fundamental role in distinguishing different ground states.
In Section 5.2, we derive our main result concerning these models. We consider local translation-invariant perturbations to the Hamiltonian of such a model, and show that the effective Hamiltonian is a linear combination of string-operators, i.e., up to an irrelevant global energy shift. The coefficients {f a } a are determined by the perturbation. They can be expressed in terms of a certain sum of diagrams, as we explain below.
While not essential for our argument, translation-invariance allows us to simplify the parameter dependence when expressing the coefficients f a and may also be important for avoiding the proliferation of small gaps. We emphasize that the effective Hamiltonian has the form (15) independently of the choice of perturbation. The operators {F a } a are mutually commuting, and thus have a distinguished simultaneous eigenbasis (we give explicit expressions for the latter in Section 5.1.6). The effective Hamiltonian (15) is therefore diagonal in a fixed basis irrespective of the considered perturbation. Together with the general reasoning for Conjecture 1, this suggests that Hamiltonian interpolation can only prepare a discrete family of different ground states in these anyonic systems.
In Section 6, we consider two-dimensional topologically ordered systems and find effective Hamiltonians analogous to (15). We will also show numerically that Hamiltonian interpolation indeed prepares corresponding ground states.

Background on anyon chains
The models we consider here describe effective degrees of freedom of a topologically ordered system. Concretely, we consider one-dimensional chains with periodic boundary conditions, where anyonic excitations may be created/destroyed on L sites, and may hop between neighboring sites. Topologically (that is, the language of topological quantum field theory), the system can be thought of as a torus with L punctures aligned along one fundamental cycle. Physically, this means that excitations are confined to move exclusively along this cycle (we will consider more general models in section 6). A well-known example of such a model is the Fibonacci golden chain [FTL + 07]. Variational methods for their study were developed in [PCB + 10, KB10], which also provide a detailed introduction to the necessary formalism. In this section, we establish notation for anyon models and review minimal background to make the rest of the paper self-contained.

Algebraic data of anyon models: modular tensor categories
Let us briefly describe the algebraic data defining an anyon model. The underlying mathematical object is a tensor category. This specifies among other things: (i) A finite set of particle labels A = {1, a, . . .} together with an involution a →ā (called particle-anti-particle exchange/charge conjugation). There is a distinguished particle 1 = 1 called the trivial or vacuum particle.
(ii) A collection of integers N c ab indexed by particle labels, specifying the so-called fusion multiplicities (as well as the fusion rules). For simplicity, we will only consider the multiplicityfree case, where N c ab ∈ {0, 1} (this captures many models of interest). In this case, we will write N c ab = δ abc . (iii) A 6-index tensor F : A 6 → C (indexed by particle labels) F abe cdf which is unitary with respect to the rightmost two indices (e, f ) and can be interpreted as a change of basis for fusion trees.
(iv) A positive scalar d a for every particle label a, called the quantum dimension.
(v) A unitary, symmetric matrix S ij indexed by particle labels such that Sī j = S ij .
(vi) A topological phase e iθ j , θ j ∈ R, associated with each particle j. We usually collect these into a diagonal matrix T = diag({e iθ j } j ); the latter describes the action of a twist in the mapping class group representation associated with the torus (see Section 6.2).
A list of the algebraic equations satisfied by these objects can be found e.g., in [LW05] (also see [NSS + 08, LW05, Kit06, Wan10] for more details). Explicit examples of such tensor categories can also be found in [LW05], some of which we discuss in Section 6.3.2.
Here we mention just a few which will be important in what follows: the fusion rules δ ijk are symmetric under permutations of (i, j, k). They satisfy m δ ijm δ mk¯ = m δ jkm δ im¯ which expresses the fact that fusion (as explained below) is associative, as well as Some of the entries of the tensor F are determined by the fusion rules and the quantum dimensions, that is, Another important property is the Verlinde formula which is often summarized by stating that S "diagonalizes the fusion rules".

The Hilbert space
The Hilbert space of a one-dimensional periodic chain of L anyons is the space associated by a TQFT to a torus with punctures. It has the form where the indices a j , b k are particle labels, V ab c are the associated finite-dimensional fusion spaces and we identify b 0 = b L . The latter have dimension dim V ab c = N c ab . Again, we will focus on the multiplicity-free case where N c ab = δ abc ∈ {0, 1}. In this case, we can give an orthonormal where a = (a 1 , . . . , a L ) and b = (b 1 , . . . , b L ) have to satisfy the fusion rules at each vertex, i.e., dim V The prefactor in the definition of the state (19) involves the quantum dimensions of the particles, and is chosen in such a way that {| a, b } is an orthonormal basis with respect to the inner product defined in terms of the isotopy-invariant calculus of diagrams: the adjoint of | a, b is represented as .

Inner products and diagramatic reduction rules
Inner products are evaluated by composing diagrams and then reducing, i.e., where [·] vac is the coefficient of the empty diagram when reducing. Reduction is defined in terms of certain local moves. These include (i) reversal of arrows (together particle-antiparticle involution a →ā) (ii) (arbitrary) insertions/removals of lines labeled by the trivial particle 1. Since1 = 1, such lines are not directed, and will often be represented by dotted lines or omitted altogether, which leads to a formal linear combination of diagrams where subgraphs are replaced locally by the figure on the rhs.
(iv) removal of "bubbles" by the substitution rule These reduction moves can be applied iteratively in arbitrary order to yield superpositions of diagrams. An important example of this computation is the following: The series of steps first makes use of an F -move (21), followed by Eq. (17) as well as (22). Together with property (16) and evaluation of the inner product (20), this particular calculation shows that the flux-eigenstates (27) are mutually orthogonal. We refer to [LW05] for more details.

Local operators
Operators are also defined by diagrams, and are applied to vectors/multiplied by stacking (attaching) diagrams on top of the latter. Expressions vanish unless all attachment points have identical direction and labels. Here we concentrate on 1-and 2-local operators, although the generalization is straightforward (see [KB10,Bon09]). A single-site operatorĤ is determined by coefficients { a } a and represented at It acts diagonally in the fusion tree basis, i.e., writing H j for the operatorĤ applied to site j, we have A two-site operatorV acting on two neighboring sites is determined by a tensor {α rs ef g } r,s,e,f,g (where the labels have to satisfy appropriate fusion rules) via the linear combinations of diagrams 21 When applied to sites j and j + 1 it acts as where the rhs. specifies a vector in H in terms of the reduction rules. It will be convenient in the following to distinguish between linear combinations of the form (24) and operators which are scalar multiplies of a single diagram (i.e., with only one non-zero coefficient α rs ef g ). We call the latter kind of two-site operator elementary.
We can classify the terms appearing in (24) according to the different physical processes they represent: in particular, we have pair creation-and annihilation operators simultaneous annihilation-and creation operatorŝ left-and right-moving 'propagation' terms as well as more general fusion operators such as e.g.,āb (We are intentionally writing down a linear combination here.) Note that a general operator of the form (24) also involves braiding processes since can be resolved to diagrams of the form using the R-matrix (another object specified by the tensor category). We will consider composite processes composed of such two-local operators in Section 5.1.7.

Ground states of anyonic chains
We will consider translation-invariant Hamiltonians H 0 = jĤ j with local terms of the form with a > 0 for a = 1 and 1 = 0 . 22 Such a Hamiltonian H 0 corresponds to an on-site potential for anyonic excitations, where a particle of type a has associated energy a independently of the site j. We denote the projection onto the ground space of this Hamiltonian by P 0 . This is the space where 1 = (1, . . . , 1) and b · 1 = (b, . . . , b). In other words, the ground space of H 0 is degenerate, with degeneracy equal to the number of particle labels. It will be convenient to use the basis {|b } b of the ground space consisting of the 'flux' eigenstates In addition, we can define a dual basis {|b } b of the ground space using the S-matrix. The two bases are related by for all particle labels a, b.
As we discuss in Section 6.3.2, in the case of two-dimensional systems, the dual basis (28) is simply the basis of flux eigenstates with respect to a 'conjugate' cycle. While this interpretation does not directly apply in this 1-dimensional context, the basis {|a } a is nevertheless welldefined and important (see Eq. (30)).

Non-local string-operators
In the following, certain non-local operators, so-called string-operators, will play a special role. Strictly speaking, these are only defined on the subspace (26). However, we will see in Section 5.2 that they arise naturally from certain non-local operators.
The string-operators {F a } a are indexed by particle labels a. In terms of the basis (27) of the ground space P 0 H of H 0 , the action of F a is given in terms of the fusion rules as Proof. We first expand P 0 into its span and F b according to eq. (29), followed by an expansion of N d bc through the Verlinde formula (18). Finally, we use the unitarity and symmetry of S to transform bra and ket factors into the dual basis given by Eq. (28)

Products of local operators and their logical action
Operators preserving the ground space P 0 H (cf. (27)) are called logical operators. As discussed in Section 5.1.6, string-operators {F a } are an example of such logical operators. Clearly, because they can simultaneously be diagonalized (cf. (30)), they do not generate the full algebra of logical operators. Nevertheless, they span the set of logical operators that are generated by geometrically local physical processes preserving the space P 0 H. That is, if O = j k V j,k is a linear combinations of products of local operators V j,k , then its restriction to the ground space is of the form i.e., it is a linear combination of string operators (with some coefficients o a : This processes has trivial action on the ground space: it is entirely local. It has action P 0 O 1 P 0 = d a P 0 , where the proportionality constant d a results from Eq. (22).
This process creates a particles anti-particle pair (a,ā) and further separates these particles. Since the operator maps ground states to excited states, This process involves the creation of a pair of particles (a,ā), with subsequent propagation and annihilation. Its logical action is P 0 O 2 P 0 = F a is given by the string-operator F a , by a computation similar to that of (23).

Perturbation theory for an effective anyon model
In this section, we consider a 1-dimensional translation-invariant system of anyons described by the Hamiltonian H 0 introduced in (25). We further consider a translation-invariant two-local perturbation V = jV j,j+1 with local termsV j,j+1 of the form (24) given bŷ where V R collects all other two-anyon processes (it will turn out that in lowest order perturbation theory, only creation and propagation are relevant). The choice of complex conjugate pairs of parameters ensures that the perturbation is self-adjoint. We may think of γ a as the 'creation ampitude', τ a as the 'propagation amplitude', and a as the energy of particle a. We now compute the form of the effective Schrieffer-Wolff-Hamiltonian. Our main result is the following: a a a a aV C j,j+1 (a) corresponds to a process where a particle pair (a,ā) is created, there is some propagation, and the particles fuse subsequently. This has trivial action on the ground space, i.e., ,2 corresponds to a process where a pair (a,ā) of particles is created, and they propagate all the way around the chain before annihilating. Its action on the ground space is given by the string-operator Lemma 5.1 (Effective Hamiltonians for 1-dimensional anyon chains). Consider H 0 + V , with the perturbation V as described. Let P 0 be the projection onto the ground space of H 0 . Then the L-th order effective Hamiltonian has the form for some constant c ∈ R, and some function f L which is independent of the particle label a and is a homogeneous polynomial of degree L in γ a and τ a .
Clearly, the form Eq. (33) of the effective Hamiltonian is consistent with the topological superselection rule (31). However, Eq. (33) provides additional information: for example, the coefficient of the string-operator F a only depends on the energy a of anyon a, as well as its creation/annihilation (γ a respectively γ a ) and propagation (τ a ) amplitudes. There is no dependence on particles distinct from a (and corresponding braiding processes). Such terms only enter in higher orders of the perturbative series. This can be thought of as a rigorous derivation of the tunneling amplitude for a particle in the weak perturbation limit. We note that due to f L being homogeneous of degree L, the dominant tunneling process will be highly sensitive to the perturbation strengths associated to different anyon labels a for large system sizes L. In the absence of a symmetry or fine tuning, it should be possible to order the terms f L ( a , γ a , τ a ) by absolute value, with different orders of magnitude being expected for each term (see Section 6.1 for further discussion).
Proof. It is easy to check that the conditions of Theorem 3.2 are satisfied with L equal to the length of the chain. Indeed, (L − 1)-local terms have trivial action on the ground space as discussed in Section 5.1.7. It thus suffices to consider expressions of the form Inserting the definition (32) of V , and diagrammatically expanding each term as in Section 5.1.4, we are left with a linear combination of terms of the form where V α j is a local operator given by an elementary (two-anyon) diagram (not a linear combination). Since such operators V α j map eigenstates of H 0 to eigenstates, and the energies of excited states reached from the ground space by applying such operators is independent of the ground state considered, each operator G merely adds a scalar, i.e., we have for some constant θ depending on the perturbations {V α j }. But the rhs. of this equation is a product of local operators as considered in Section 5.1.7. According to the expression (31), this is a linear combination of string-operators, i.e., Furthermore, since each V α j is an elementary two-local operator, and we consider only products of length L, the only terms P 0 V α 1 V α 2 V α 3 · · · V α L P 0 that have non-trivial action on the ground space are those associated with processes where a single particle (say of type a) winds around the whole chain. We will call such a process topologically non-trivial. Its action on the ground space is given by a single string-operator F a .
In summary (rearranging the sum), we conclude that the L-th order effective Hamiltonian has the form (33), where the coefficient f L ( a , γ a , τ a ) has the form and where the sum is over the set of all length-L-topologically non-trivial processes (consisting of elementary terms) involving particle a. The coefficient Furthermore, ν(V α 1 , . . . , V α L ) can only be non-zero when all L operators V α j are either pair creation/anihilation or hopping terms involving the particle a. This implies the claim.

2D topological quantum field theories
As discussed in Section 4, adding a local perturbation to a Majorana chain leads to an effective Hamiltonian given by the parity (string)-operator. Similarly, in the case of a general anyon chain (discussed in Section 5), the effective Hamiltonian is a linear combination of stringoperators F a , associated with different particle labels a. Here we generalize these considerations to arbitrary systems described by a 2-dimensional topological quantum field theory (TQFT) and subsequently specialize to microscopic models, including the toric code and the Levin-Wen string-net models [LW05].
Briefly, a TQFT associates a "ground space" H Σ to a two-dimensional surface Σ -this is e.g., the ground space of a microscopic model of spins embedded in Σ with geometrically local interactions given by some Hamiltonian H 0 (see Section 6.3). In other words, H Σ ⊂ H phys,Σ is generally a subspace of a certain space H phys,Σ of physical degrees of freedom embedded in Σ. The system has localized excitations (anyons) with (generally) non-abelian exchange statistics. In particular, there are well-defined physical processes involving creation, propagation, braiding and annihilation of anyons, with associated operators as in the case of 1-dimensional anyon chains (see Section 5). Contrary to the latter, however, the particles are not constrained to move along a 1-dimensional chain only, but may move arbitrarily on the surface Σ. Nevertheless, the description of these processes is analogous to the case of spin chains, except for the addition of an extra spatial dimension. For example, this means that local operators acting on a region R ⊂ Σ are now represented by a linear combination of string-nets (directed trivalent graphs with labels satisfying the fusion rules) embedded in R×[0, 1]. We refer to e.g., [FKLW03] for more examples of this representation.
As before, there are distinguished ground-space-to-ground-space (or "vacuum-to-vacuum") processes which play a fundamental role. These are processes where a particle-anti-particle pair (a,ā) is created, and the particles fuse after some propagation (tunneling), i.e., after tracing out a closed loop C on Σ. Non-trivial logical operators must necessarily include topologically non-trivial loops C on Σ in their support (the spatial region in which they are physically realized). In particular, for any such loop C, there is a collection {F a (C)} a of string-operators associated with different particle labels. More precisely, a loop is a map C : [0, 1] → Σ satisfying C(0) = C(1). Reversing direction of the loop gives a new loopC(t) := C(1 − t), and this is equivalent to interchanging particle-and antiparticle labels: we have the identity F a (C) = Fā(C). In Section 6.2, we state some general properties of the string-operators {F a (C)} a , and, in particular, explain how to express them in suitable bases of the ground space.

Perturbation theory for Hamiltonians corresponding to a TQFT
In general, the anyon model associated with a TQFT is emergent from a microscopic spin Hamiltonian H 0 . The anyon notion of site, as discussed in Section 5, does not necessarily coincide with the spin notion of site associated with the microscopic spin model. Nevertheless, the following statements are true: (i) any non-trivial logical operator must include at least one non-contractible loop in its support.
(ii) given a perturbation V consisting of geometrically local operators, there exists some minimum integer L such that H 0 , V satisfy the topologically ordered condition with parameter L.
In general, the value of L will depend on the length of the shortest non-contractible loop(s), and the resulting effective Hamiltonian will be of the form where the dependence on H 0 and the coefficients in V has been left implicit. The sum is over all non-trivial loops C of length L (where length is defined in terms of the spin model), as well as all particle labels a.
Computing the coefficients {f L (a, C)} may be challenging in general. Here we discuss a special case, where anyon processes associated with a single particle a (respectively its antiparticleā) are dominant (compared to processes involving other particles). That is, let us assume that we have a translation-invariant perturbation V of the form where the sum is over all pairs (j, j ) of nearest-neighbor (anyonic) sites, andV (1) j,j =V (1) and V (2) j,j =V (2) are both 1-and 2-local operators on the same anyon site lattice -this is a straightforward generalization of anyon chains to 2D. Our specialization consists in the assumption that all local creation, propagation and annihilation processes constituting the operatorV (1) j,j =V (1) only correspond to a single anyon type a (andā), and that these processes are dominant in the sense that the remaining terms satisfy ηV (2) V (1) . In the limit η → 0, perturbation theory in this model only involves the particles (a,ā).

28
Assuming that the shortest non-contractible loops have length L in this anyonic lattice, we claim that where G (L) eff is an effective Hamiltonian with the same form as H (L) eff ( ), but only contains string operators F b (C) with b = a. The reason is that in order to generate a string operator F b (C) in L steps (i.e., at L-th order in perturbation theory), we need to apply local operators corresponding to anyon b L times, as discussed in Lemma 5.1. Such local operators can only be found in ηV 2 , therefore we obtain the coefficient eff . Thus if we fix the system size and slowly increase η from 0, the (relative) change of the total effective Hamiltonian is exponentially small with respect to L. This implies that the ground state of the effective Hamiltonian is stable when η is in a neighbourhood of 0. We will see in Section 7 that the final states of Hamiltonian interpolation are indeed stable in some regions of initial Hamiltonians. The above discussion can be viewed as a partial explanation 6 for this phenomenon.

String-operators, flux bases and the mapping class group
In the following, we explain how to compute effective Hamiltonians of the form (35) in the case where the perturbation is isotropic, resulting in identical coefficients f L (a, C) = f L (a, C ) for all loops C of identical length. This will be guaranteed by symmetries. We give explicit examples in Section 7.
For this purpose, we need a more detailed description of the action of string-operators on the ground space. Consider a fixed (directed) loop C : [0, 1] → Σ embedded in the surface Σ. The process of creating a particle-anti-particle pair (a,ā), then propagating a along C, and subsequently fusing withā defines an operator F a (C) which preserves the ground space H Σ . The family of operators {F a (C)} a is mutually commuting and defines a representation of the Verlinde algebra. It is sometimes convenient to consider the associated (images of the) idempotents, which are explicitly given by (as a consequence of the Verlinde formula (18)) The operators P a (C) are mutually orthogonal projections P a (C)P b (C) = δ ab P a (C). The inverse relationship (using the unitarity of S) reads and is the generalization of (30): indeed, specializing to the case where Σ is the torus (this will be our main example of interest), and C is a fundamental loop, the operators P a (C) are rank-one projections (when restricted to the ground space), and determine (up to phases) an orthonormal basis of B C = {|a C } a of H Σ by P a (C) = |a C a C |. In physics language, the state |a C has "flux a" through the loop C. (More generally, one may define "fusion-tree" basis for higher-genus surfaces Σ by considering certain collections of loops and the associated idempotents, see e.g., [KKR10]. However, we will focus on the torus for simplicity.) Consider now a pair of distinct loops C and C . Both families {F a (C)} a and {F a (C )} a of operators act on the ground space, and it is natural to ask how they are related. There is a simple relationship between these operators if C = ϑ(C) is the image of C under an element ϑ : Σ → Σ of the mapping class group MCG Σ of Σ (i.e., the group of orientation-preserving diffeomorphisms of the surface): The TQFT defines a projective unitary representation V : MCG Σ → U(H Σ ) of this group on H Σ , and we have In general, while the topology of the manifold is invariant under the mapping class group, the specific lattice realization may not be. For this reason, if we desire to lift the representation V to the full Hilbert space H Σ ⊃ H phys,Σ , such that the resulting projective unitary representation preserves the microscopic Hamiltonian H 0 under conjugation, we may need to restrict to a finite subgroup of the mapping class group MCG Σ . If the lattice has sufficient symmetry, such as for translation-invariant square or rhombic lattices, one may exploit these symmetries to make further conclusions about the resulting effective Hamiltonians.

String-operators and the mapping class group for the torus
For the torus, the mapping class group MCG Σ is the group SL(2, Z). To specify how a group element maps the torus to itself, it is convenient to parametrize the latter as follows: we fix complex numbers (e 1 , e 2 ) and identify points z in the complex plane according to z ≡ z + n 1 e 1 + n 2 e 2 for n 1 , n 2 ∈ Z .
In other words, (e 1 , e 2 ) defines a lattice in C, whose unit cell is the torus (with opposite sides identified). A group element A = a b c d ∈ SL(2, Z) then defines parameters (e 1 , e 2 ) by e 1 = ae 1 + be 2 e 2 = ce 1 + de 2 , which a priori appear to be associated with a new torus. However, the constraint that A ∈ SL(2, Z) ensures that (e 1 , e 2 ) and (e 1 , e 2 ) both define the same lattice, and this therefore defines a map from the torus to itself: The action of A is given by αe 1 + βe 2 → αe 1 + βe 2 for α, β ∈ R, i.e., it is simply a linear map determined by A.
The group SL(2, Z) = t, s is generated by the two elements t = Dehn twist 1 1 0 1 and π/2 rotation s = 0 1 −1 0 which are equivalent to the Möbius transformations τ → τ +1 and τ → −1/τ . Clearly, t fixes e 1 and hence the loop C : t → C(t) = te 1 , t ∈ [0, 1] on the torus (this loop is one of the fundamental cycles). The matrices representing the unitaries V (t) and V (s) in the basis B C = {|a C } a of H Σ (where |a C is an eigenstate of P a (C) = |a C a C |) are denoted T and S, respectively. These e 1 e 2 C 2 C 1 Figure 2: Minimal loops on the square torus matrices are given by the modular tensor category: T is a diagonal matrix with T aa = e iθa (where θ a is the topological phase of particle a), whereas S is the usual S-matrix. This defines the mapping class group representation on the Hilbert space H Σ associated with the torus Σ.
In the following, we compute explicit relationships between string-operators of minimal length. We consider two cases: a square torus and a rhombic torus. This allows us to express terms such as those appearing in Eq. (34) in a fixed basis.
Square torus. Here we have e 1 = 1 and e 2 = i .
There are (up to translations) two loops of minimal length, which may be traversed in either of two directions namely for t ∈ [0, 1], see Fig. 2. Since se 1 = −e 2 and se 2 = e 1 , we conclude that C 2 (t) = s(C 1 (t)) C 1 (t) = s 2 (C 1 (t)) C 2 (t) = s 3 (C 1 (t)) C 1 (t) = s 4 (C 1 (t)) In particular, expressed in the basis B C 1 , we have Thus, when the lattice and Hamiltonian H 0 obey a π/2 rotation symmetry, the effective perturbation Hamiltonian will be proportional to (38). This is the case for the toric code on a square lattice.

Minimal loops of interest are shown in
e 1 e 2 C 1 C 2 C 3 Figure 3: Minimal loops on the rhombic torus for t ∈ [0, 1]. Observe that these can be related by a π/3 rotation u (if we use the periodicity of the lattice), i.e., Since such a rotation u maps e 1 , e 2 to e 1 = e 2 e 2 = e 2 − e 1 , it is realized by the element u = 0 1 −1 1 ∈ SL(2, Z), which decomposes into the generators (37) as u = ts 3 ts. We conclude that, expressed in the basis B C 1 , we have Again, if the lattice and Hamiltonian H 0 are invariant under a π/3 rotation, we may conclude that the effective perturbation Hamiltonian will have the form (39). This is the case for the Levin-Wen model on a honeycomb lattice embedded in a rhombic torus (see also Section 7.2).

Microscopic models
The purpose of this section is two-fold: First, we briefly review the construction of the microscopic models we use in our numerical experiments in Section 7: these include the toric code (see Section 6.3.1) as well as the doubled semion and the doubled Fibonacci model, both instantiations of the Levin-Wen construction (see Section 6.3.2). Second, we define single-qubit operators in these models and discuss their action on quasi-particle excitations (i.e., anyons). This translation of local terms in the microscopic spin Hamiltonian into operators in the effective anyon models is necessary to apply the perturbative arguments presented in Section 6.1. We will use these local terms to define translation-invariant perturbations (respectively trivial initial Hamiltonians) in Section 7).

The toric code
Kitaev's toric code [Kit03] is arguably the simplest exactly solvable model which supports anyons. It can be defined on a variety of lattices, including square and honeycomb lattices.
Here we will introduce the Hamiltonian corresponding to honeycomb lattice. On each edge of the lattice resides a qubit. The Hamiltonian consists of two parts and takes the form where B p = X ⊗6 is the tensor product of Pauli-X operators on the six edges of the plaquette p, and A v = Z ⊗3 is the tensor product of Pauli-Z operators on the three edges connected to the vertex v.
Note that in terms of its anyonic content, the toric code is described by the double of Z 2 ; hence a model with the same type of topological order could be obtained following the prescription given by Levin and Wen (see Section 6.3.2). Here we are not following this route, but instead exploit that this has the structure of a quantum double (see [Kit03]). The resulting construction, given by (40), results in a simpler plaquette term B p as opposed to the Levin-Wen construction.
The anyonic excitations supported by the toric code are labeled by {1, e, m, }. The e anyon or electric excitation corresponds to vertex term excitations. The m anyon or magnetic excitations correspond to plaquete term excitations. Finally, the anyon corresponds to an excitation on both plaquete and vertex and has the exchange statistics of a fermion. We can write down the string operators F a (C) for a closed loop C on the lattice explicitly (see [Kit03]). Without loss of generality, we can set F e (C) = P 0 i∈C X i P 0 and F m (C) = P 0 i∈D Z i P 0 , where D is a closed loop on the dual lattice corresponding to C. Finally, the operator F (C) = F e (C) × F m (C) can be written as a product of F e (C) and F m (C), since e and m always fuse to . With respect to the ordering (1, e, m, ) of the anyons, the S-and T -matrices described in Section 5.1.1 are given by for the toric code.
Local spin operators. A natural basis of (Hermitian) operators on a single qubit is given by the Pauli operators. For the toric code, each of these operators has a natural interpretation in terms of the underlying anyon model. Consider for example a single-qubit Z-operator. The "anyonic lattice" associated with manyons is the dual lattice (i.e., these anyons 'live' on plaquettes), and a single-qubit Z-operator acts by either creating or annihilating a (m,m) = (m, m) on the neighboring plaquettes, or propagating an existing m from one plaquette to the other. That is, in the terminology of Section 5.1.4, a Z-operator acts as a local term in the effective anyon model. An analogous identity holds for X, which is associated with eanyons: the latter live on vertices of the spin lattice. Finally, Y -operators act on -anyons in the same manner; these anyons live on 'supersites', consisting of a plaquette and and an adjacent vertex.

Short introduction to the Levin-Wen model
Levin and Wen [LW05] define a family of frustration-free commuting Hamiltonian with topologically ordered ground space and localized anyonic excitations. Their construction is based on interpreting the state of spins residing on the edges of a trivalent lattice (such as a honeycomb lattice) as configurations of string-nets.
To specify a string-net model, we need algebraic data associated with an anyon model as described in Section 5.1.1. This specifies, in particular, a set of anyon labels F = {a i }, associated fusion rules, as well as S-and F -matrices. The Levin-Wen model then associates a qudit to each edge of the lattice, where the local dimension of each spin corresponds to the number of anyon labels in F. One chooses an orthonormal basis {|a } a∈F ⊂ C |F | indexed by anyon labels; in the following, we usually simply write a instead of |a to specify a state of a spin in the microscopic model. The Levin-Wen spin Hamiltonian can be divided into two parts, where each B p is a projector acting on the 12 edges around a plaquette p, and each A v is a projector acting on the 3 edges around a vertex v. In particular, we can construct the spin Hamiltonian for the doubled semion and the doubled Fibonacci models in this way by choosing different initial data.
As long as all the particles in the underlying model F are their own antiparticles (i.e., the involution a →ā is the identity), it is not necessary to assign an orientation to each edge of the lattice. This affords us an important simplification, which is justified for the models under consideration: these only have a single non-trivial anyon label, which is itself its own antiparticle (recall that the trivial label satisfies1 = 1). With this simplification, which we will use throughout the remainder of this paper, the vertex operator A v can be written as where d s is the quantum dimension of the anyon label s, and D = j d 2 j is the total quantum dimension.
Having specified the spin Hamiltonian, we stress that the anyon labels F used in this construction should not be confused with the anyon labels D(F) describing the local excitations in the resulting Hamiltonian (43). The latter can be described as 'pairs' of anyons from F, i.e., D(F) = {(a i , a j )} a i ,a j ∈F . Their fusion, twist and braiding properties are described by the double of the original theory. The S D(F ) -and F D(F ) -matrices of D(F) can be obtained from the S-and T -matrix associated with F (see [LW05]). String operators F a i ,a j (C) acting on the spin lattice have also been explicitly constructed in [LW05] Below, we present some of the specifics of two models constructed in this way: the doubled semion and doubled Fibonacci model. In addition to Kitaev's toric codes D(Z 2 ), these are the only models defined on two labels (i.e., with microscopic qubit degrees of freedom).

The doubled semion model
The underlying string-net model of the doubled semion model only consists of one non-trivial label s and the trivial label 1. To specify the spin Hamiltonian, we have d s = 1, and δ abc = 1 if and only if an even number of a, b, c are s. The F -matrix is given by F ss1 ss1 = −1 and otherwise F abc def is 0 or 1 depending on whether (a, b, c, d, e, f ) is a legal configuration (see [Kit06] for more detailed explanation). As we explained above, to construct a spin Hamiltonian, we put a qubit on each edge of the lattice with orthonormal basis |1 , |s . The spin Hamiltonian obtained this way is similar to the toric code and it also supports Abelian anyons. The excitations of the spin model can be labeled by D(F) = {(1, 1), (1, s), (s, 1), (s, s)}, which is the quantum double of F = {1, s}. With respect to the given ordering of anyons, the S-and T -matrices of these excitations are given by Local operators. Identifying |1 with the standard basis state |0 and |s with |1 , we can again use Pauli operators to parametrize single-spin Hamiltonian terms. Here we will discuss the effect of single qubit operators X and Z on the ground states of the resulting topologically ordered Hamiltonian. The goal is to interpret single spin operators in terms of effective anyon creation, annihilation and hopping operators.
When Z-operator is applied to an edge of the system in a ground state, only the neighboring plaquete projectors B p will become excited. More specifically, a pair of (s, s) anyons are created if none were present. Since (s, s) is an abelian anyon, in fact a boson, and is the anti-particle of itself, a Z operator could also move an (s, s) anyon or annihilate two such particles if they are already present. Thus we conclude that single-qubit Z-operators have a similar action as in the toric code (cf. (42)), with s playing the role of the anyon m.
When an X operator is applied on edge of the system in a ground state, it excites the two neighboring vertex terms A v (in the sense that the state is no longer a +1-eigenstate any longer). Since the plaquete terms B p are only defined within the subspace stabilized by A v , the four plaquette terms B p terms around the edge also become excited. It is unclear how to provide a full interpretation of X operators in terms of an effective anyon language. In order to provide this, a full interpretation of the spin Hilbert space and its operators in the effective anyonic language is required; such a description is currently not known.
In summary, this situation is quite different from the case of the toric code, where X and Z are dual to each other.

The doubled Fibonacci
Again, the underlying string-net model of doubled Fibonacci contains only one non-trivial label τ , with quantum dimension d τ = ϕ, where ϕ = 1+ √ 5 2 . The fusion rules are given by δ abc = 0 if only one of the a, b, c is τ , and otherwise δ abc = 1. Non-trivial values of F are and otherwise F abc def is either 0 or 1 depending on whether (a, b, c, d, e, f ) is a legal configuration. Many aspects of the doubled Fibonacci spin Hamiltonian are similar to the doubled semion model: • There is one qubit on each edge, with orthonormal basis states associated with the anyon labels F = {1, τ }. With respect to the given ordering of anyons, the S-and T -matrices are given by A substantial difference to the doubled semion model is that the non-trivial anyons supported by the model are non-abelian. One manifestation of this fact we encounter concerns the (τ , τ )-anyon: • While (τ , τ ) is its own anti-particle, it is not an abelian particle so in general two (τ , τ ) particles will not necessarily annihilate with each other. In other words, the dimension of the subspace carrying two localized (τ , τ ) charges is larger than the dimension of the charge-free subspace.
Neither of these properties holds for the (s, s)-anyon in the case of the doubled semion model.
Local operators. Similarly, as before, we identify |1 with the standard basis state |0 and |τ with |1 , enabling us to express single-qubit operators in terms of the standard Pauli operators. Again, we want to consider the effect of single qubit operators in terms of anyons. This is generally rather tricky, but for single-qubit Z-operators, we can obtain partial information from an analysis presented in appendix B: Let |ψ be a ground state. Then Z|ψ = 1

Numerics
In this section, we present results obtained by numerically simulating Hamiltonian interpolation for small systems. Specifically, we consider three topologically ordered systems on the 12-qubit honeycomb lattice of Fig 4: the toric code, the doubled semion and the doubled Fibonacci Levin-Wen models. That is, the target Hamiltonian H top is given either by (40) (with stabilizer plaquette-and vertex-operators A v and B p ) in the toric code case, and expression (43) specified in Section 6.3.2 (with projection operators A v and B p ) for the doubled semion and the doubled Fibonacci case. As initial Hamiltonian H triv , we choose certain translation-invariant Hamiltonians consisting of single-qubit Pauli-X, Pauli-Y and Pauli-Z operators (see Sections 6.3.3 and 6.3.4 for their definition and a discussion of the effect of these operators in the two Levin-Wen models.) For concreteness and ease of visualization, we will consider the following families of such Hamiltonians: the one-parameter family where θ ∈ [0, 2π[, and two two-parameter families of the form where (a, b) ∈ R 2 belongs to the unit disc, a 2 + b 2 ≤ 1. (In some instances, we will permute the roles of X, Y and Z, and use an additional superscript to indicate this.) For different parameter choices θ respectively (a, b), we study Hamiltonian interpolation (i.e., the evolution (2)) along the linear interpolation path H(t) (cf. (1)) with a total evolution time T . In order to numerically simulate the evolution under the time-dependent Schrödinger equation, we perform a time-dependent Trotter expansion using the approximation Unless otherwise specified, the time discretization is taken to be ∆t = 0.1.

Quantities of interest and summary of observations
Recall that our initial state Ψ(0) = ϕ ⊗12 is the unique 12-qubit ground state of the chosen trivial Hamiltonian H triv . We are interested in the states Ψ(t) along the evolution, and, in particular, the final state Ψ(T ) for a total evolution time T . For notational convenience, we will write Ψ θ (t), respectively Ψ ± a,b (t) to indicate which of the initial Hamiltonians H triv is considered (cf. (46) and (47)). We consider the following two aspects: (Non)-adiabaticity: We investigate whether the state Ψ(t) follows the instantanenous ground space along the evolution (2). We quantify this using the adiabaticity error, which we define (for a fixed total evolution time T , which we suppress in the notation) as where P 0 (t) is the projection onto the ground space of H(t) (note that except for t = T , where P 0 (T ) projects onto the degenerate ground space of H top , this is generally a rankone projection). The function t → adia (t) quantifies the overlap with the instantaneous ground state of H(t) along the Hamiltonian interpolation t → H(t), and hence directly reflects adiabaticity.
Ultimately, we are interested in whether the evolution reaches a ground state of H top . This is measured by the expression adia (T ), which quantifies the deviation of the final state Ψ(T ) from the ground space of H top . Clearly, the quantity adia (T ) depends on the choice of initial Hamiltonian H triv (i.e., the parameters θ respectively (a, b)) and the total evolution time T . For sufficiently large choices of the latter, we expect the adiabaticity assumption underlying Conjecture 1 to be satisfied, and this is directly quantifiable by means of the adiabaticity error. We will also discuss situations where, as discussed in Observation 2.3, symmetries prevent reaching the ground space of H top as reflected in a value of adia (T ) close to 1.
Logical state: assuming the ground space of H top is reached (as quantified by adia (T )), we will identify the logical state Ψ(T ) and investigate its stability under perturbations of the the initial Hamiltonian H triv (i.e., changes of the parameters θ respectively (a, b)). For this purpose, we employ the following measures: • We argue (see Section 7.2) that symmetries constrain the projection of the resulting state Ψ(T ) onto the ground space of H top to a two-dimensional subspace (see Section 7.2). For the toric code, the state is then fully determined by the expectation values X Ψ(T ) , Z Ψ(T ) of two logical operatorsX andZ. To investigate stability properties of the prepared state, we can therefore consider ( X Ψ(T ) , Z Ψ(T ) ) as a function of parameters of the initial Hamiltonian.
• for the Levin-Wen models, we proceed as follows: we pick a suitable reference state |ψ R ∈ (C 2 ) ⊗12 in the ground space of H top , and then study how the overlap | Ψ ± a,b (T )|ψ R | 2 changes as the parameters (a, b) of the initial Hamiltonian are varied. In particular, if we fix a pair (a 0 , b 0 ) and choose |ψ R as the normalized projection of the state |Ψ ± a 0 ,b 0 (T ) onto the ground space of H top , this allows us to study the stability of the prepared state |Ψ ± a,b (T ) as a function of the Hamiltonian parameters (a, b) in the neighborhood of (a 0 , b 0 ).
According to the reasoning in Section 3.2 (see Conjecture 1), the specific target state |ψ ref a 0 ,b 0 chosen in this way should correspond to the ground state of H top + H ± (a 0 , b 0 ) in the limit → 0 of infinitesimally small perturbations (or, more precisely, the corresponding effective Hamiltonian). Furthermore, according to the reasoning in Section 6.1, the family of effective Hamiltonians associated with H top + H ± (a, b) has a very specific form. This should give rise to a certain stability of the ground space as a function of the parameters (a, b). To support this reasoning, we numerically compute the (exact) ground state |ψ pert a,b of H top + H ± (a, b) for the choice = 0.001 (as a proxy for the effective Hamiltonian), and study the overlap | ψ pert a,b |ψ ref a 0 ,b 0 | 2 as a function of the parameters (a, b) in the neighborhood of (a 0 , b 0 ).
The results of our numerical experiments support the following two observations: • Hamiltonian interpolation is generically able to prepare approximate ground states of these topological models for sufficiently long total evolution times T .
• Specific final state(s) show a certain degree of stability with respect to changes in the initial Hamiltonian. The theoretical reasoning based on perturbation theory presented in Section 6 provides a partial explanation of this phenomenon.

A symmetry of the 12-qubit rhombic torus
As discussed in Section 6.3, the ground space of H top on a torus is 4-dimensional for the toric code, the doubled semion-and the Fibonacci model. In this section, we argue that adiabatic interpolation starting from a translation-invariant Hamiltonian (as considered here) yields states belonging to a two-dimensional subspace of this ground space, thus providing a simplification. Consider again the 12-qubit rhombic torus illustrated in Fig. 4. A π/3 rotation permuting the physical qubits according to (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) → (5, 7, 8, 6, 9, 12, 11, 10, 2, 4, 1, 3) defines a unitary U π/3 on C ⊗12 . Because of translation-invariance, this is a symmetry of the trivial Hamiltonian, U π/3 H triv U † π/3 = H triv , and it can easily be verified that for the models considered here, the unitary U π/3 also commutes with H top . Because of the product form of the initial state Ψ(0), it thus follows that U π/3 Ψ(t) = Ψ(t) along the whole trajectory t → Ψ(t) of adiabatic interpolation. In particular, the projection of the final state Ψ(T ) onto the ground space of H top is supported on the +1-eigenspaces space of U π/3 . As discussed in Section 6.2.1, a π/3-rotation of the rhombic torus corresponds to the modular transformation ts 3 ts. Since U π/3 realizes this transformations, its restriction to the ground space of H top can be computed from the T and S-matrices. That is, expressed in the flux bases discussed in Section 6.3, the action of U π/3 on the ground space is given by the matrix T S 3 T S, where (S, T ) are given by (41) for the toric code, as well as (44) and (45) for the doubled semion and Fibonacci models, respectively. The specific form of T S 3 T S or its eigenvectors is not particularly elucidating, but may be computed explicitly.
Importantly, the +1 eigenspace of T S 3 T S is two-dimensional for the toric code, the doubled semion and the Fibonacci models. (In the case of the toric code, it can be verified that this eigenspace is contained in the logical symmetric subspace. The latter is the subspace invariant under swapping the two logical qubits in the standard computational basis.) As a result, the projection of the state Ψ(T ) onto the ground space of H top belongs to a known two-dimensional subspace which can be explicitly computed. This means that we may characterize the resulting state in terms of a restricted reduced set of logical observables, a fact we will exploit in Section 7.3.

The toric code
As discussed in Section 6.3.1, for the toric code on the honeycomb lattice (see Fig. 4), the Hamiltonian of the model is is a tensor product of Pauli-X operators on the six edges of the plaquette p, and A v = Z ⊗3 is a tensor product of Pauli-Z operators on the three edges incident on the vertex v. We point out that the toric code on a honeycomb lattice has several differences compared to a toric code on a square lattice (which is often considered in the literature). Assuming that both lattices are defined with periodic boundary conditions, Properties (i) and (ii) imply that the usual symmetries X ↔ Z and Z ↔ −Z of the toric code on the square lattice are not present in this case. The absence of these symmetries is reflected in our simulations. Property (iii) also directly affects the final state, as can be seen by the perturbative reasoning of Section 6.1:Z-string operators appear in lower order in perturbation theory compared toX-string operators.
(Non)-adiabaticity. We first present the adiabaticity error adia (T ) for the Hamiltonian H triv (θ) given by (46) (for different values of θ) as a function of the total evolution time T . Fig. 5 illlustrates the result. It shows that for sufficiently long total evolution times T , the Hamiltonian interpolation reaches the ground space of the toric code when the initial Hamiltonian is H triv (θ = π) = − i Z i ; this is also the case for θ ∈ {π/4, π/2, 3π/4}. However, if the initial Hamiltonian is H triv (θ = 0) = i Z i , then the final state Ψ(T ) is far from the ground space of the toric code Hamiltonian H top . This phenomenon has a simple explanation along the lines of Observation 2.3. Indeed, if θ = 0, then every vertex terms A v = Z ⊗3 commutes with both H triv as well as H top (and thus all intermediate Hamiltonians H(t)). In particular, the expectation value of the vertex terms remains constant throughout the whole evolution, and this leads to an adiabaticity error adia (T ) of 1 in the case of H triv (θ = 0) = i Z i .  49)) as a function of the total evolution time T and the initial Hamiltonian chosen. For the latter, we consider the one-parameter family H triv (θ) given by (46). For θ = 0, the adiabatic evolution is not able to reach the final ground space because initially A v = −1 for every vertex operator A v = Z ⊗3 , and this quantity is conserved during the evolution. This is a feature of the honeycomb lattice because the vertex terms A v have odd weights. For other values of θ, the ground space is reached for sufficiently large total evolution times T .
• In the first case (Fig. 6a), we observe that for all initial Hamiltonians of the form H + triv (a, b) in a small neighborhood of H + triv (0, 0), the adiabaticity error adia (T ) is also large, but drops off quickly outside that neighborhood. This is consistent with the relevant level crossing(s) being avoided by introducing generic perturbations to the initial Hamiltonian.
• In contrast, almost all initial Hamiltonians in the family H − triv (a, b) (around the initial Hamiltonian H − triv (0, 0)) lead to a small adiabaticity error adia (T ) (Fig. 6b), demonstrating the stability of the adiabatic preparation.
In a similar vein, Fig. 6c illustrates the non-adiabaticity for the family of Hamiltonian The family H +X triv (a, b) (defined with a positive square root) would behave exactly the same due to the symmetry +X ↔ −X. (c) The logarithm of adiabaticity error ln adia (T ) in the neighborhood around H −X triv (0, 0) = − j X j for different Hamiltonians H −X triv (a, b). Note that the resulting figure would look identical for the Hamiltonians H +X triv (a, b) because of the −X ↔ +X symmetry. Figure 6: The adiabaticity error adia (T ) = 1 − Ψ(T )|P 0 (T )|Ψ(T ) , measuring how well the final state Ψ(T ) overlaps with the ground space of the toric code. All three figures are for a total evolution time T = 120. In Fig. 6a, we consider the family of initial Hamiltonians H + triv (a, b) in the neighborhood of H + triv (0, 0) = H triv (θ = 0) = j Z j . In contrast, Fig. 6b illustrates different choices of initial Hamiltonians H − triv (a, b) around H − triv (0, 0) = H triv (θ = π) = − j Z j . The values (a, b) ⊂ R 2 are restricted to the unit disc a 2 + b 2 ≤ 1; the center points of the two figures correspond respectively to θ = 0 and θ = π in Fig. 5. Finally, Fig. 6c gives the non-adiabaticity error for initial Hamiltonians of the form H −X triv (a, b) (as defined in Eq. (50)). 42 Logical state. For the 12-qubit rhombic toric code (Fig. 4), logical observables associated with the two encoded logical qubits can be chosen as Because of the symmetry (7.2), however, these are not independent for a state Ψ(T ) (or more precisely, its projection P 0 (T )Ψ(T )) prepared by Hamiltonian interpolation from a product state: their expectation values satisfy the identities Z 1 = Z 2 and X 1 = X 2 .
We will hence use the two (commuting) logical operators to describe the obtained logical state. In Fig. 7, we plot the expectation values ofZ andX in the final state Ψ(T ) for initial Hamiltonians of the form (cf. (47) and (50)) We again discuss the center points in more detail. It is worth noting that the single-qubit {Z i } operators correspond to the local creation, hopping and annihilation of m anyons situated on plaquettes, whereas the operators {X i } are associated with creation, hopping and annihilation of e anyons situated on vertices. In particular, this means that the initial Hamiltonians associated with the center points in the two figures each generate processes involving only either type of anyon.
• For H − triv (0, 0) = − i Z i , we know that Z = 1 during the entire evolution becauseZ commutes with the Hamiltonians H(t), and the initial ground state Ψ(0) is a +1 eigenstate ofZ. In Figs. 7a and 7b, we can see that there is a large region of initial Hamiltonians H − triv (a, b) around H − triv (0, 0) = − i Z i which lead to approximately the same final state.
• On the other hand, as shown in Figs. 7c and 7d, the stable region of Hamiltonians H −X triv (a, b) around the initial Hamiltonian H −X triv (0, 0) = − i X i is much smaller. This is due to the fact that the operatorX appears in higher order perturbation expansion compared toZ, and the evolution time T is taken to be quite long. Given sufficiently large total evolution time T , in the neighborhood of H −X triv (0, 0) = − i X i , the lower order termZ in the effective Hamiltonian will dominate the termX associated with V = − i X i . However, in both cases considered in Fig. 7, we observe that one of two specific logical states is prepared with great precision within a significant fraction of the initial Hamiltonian parameter space. (c) The quantity ln(1 − X ) for initial Hamiltonians H −X triv (a, b) in the neighborhood of H −X triv (0, 0) = − j X j . The corresponding adiabaticity error is shown in Fig. 6c.

The doubled semion model
In this section, we present our numerical results for Hamiltonian interpolation in the case of the doubled semion model (see Section 6.3.3).
(Non)-adiabaticity. We first consider the total evolution time T necessary to reach the final ground space of H top , for different initial Hamiltonians H triv . Specifically, Fig. 8 shows the adiabaticity error adia (T ) (cf. (49)) as a function of the total evolution time T for the three initial Hamiltonians H triv (θ), θ ∈ {π, π/3, 2π/3} (cf. (46)). The case of θ = 0, corresponding to the initial Hamiltonian H triv (0) = j Z j is not shown in Fig. 8 since the situation is the same as in the toric code: No overlap with the ground space of H top is achieved because the vertex-operators A v = Z ⊗3 are conserved quantities with A v = −1.
In Fig. 9a, we plot the adiabaticity error adia (T ) with initial Hamiltonian among the family of Hamiltonians H + triv (a, b) in the vicinity of H + triv (0, 0) = i Z i . Similarly, Fig. 9b provides the adiabaticity error for initial Hamiltonians H − triv (a, b) in the vicinity of H − triv (0, 0) = − i Z i .
Logical state. To explore the stability of the resulting final state, we consider the family of initial Hamiltonians H ± triv (a, b) and compute the overlap | Ψ ± a,b (T )|ψ R | 2 of the resulting final state Ψ ± a,b (T ) with a suitably chosen reference state ψ R . We choose the latter as follows: ψ R is the result of projecting the final state Ψ − 0,0 (T ) of the Hamiltonian interpolation, starting from the initial Hamiltonian H − triv (0, 0) = − i Z i onto the ground space of the doubled semion model H top and normalizing, i.e., We briefly remark that the state ψ R is uniquely determined (up to a phase) as the unique simultaneous +1-eigenvector of T S 3 T S (see Section 7.2) and the string operatorZ = Z 1 Z 2 (which is the string-operator F (s,s) (C) for the associated loop C when acting on the ground (a) The adiabaticity error adia (T ) for different Hamiltonians H + triv (a, b) in the vicinity of H + triv (0, 0) = j Z j . The adiabaticity error is maximal for the latter because of conserved quantities; however, it decays rapidly outside this center region. This situation is analogous to Fig. 6a for the toric code.  (a) The overlap | Ψ + a,b (T )|ψ R | 2 for initial Hamiltonians H + triv (a, b) around H + triv (0, 0) = i Z i . We observe that outside the center region (where the ground space of H top is not reached, see Fig. 9a), the prepared state Ψ + a,b (T ) is not too far from the reference state ψ R . Note that definition of the latter does not correspond to any Hamiltonian in this plot, but rather the centerpoint of Fig. 10b.
We plot the logarithm of this quantity because the variation is small. As illustrated, the resulting state is close to the reference state ψ R throughout most of the parameter region. Observe that, while ψ R corresponds to the center point in this figure, it still deviates from Ψ − a,b since the latter has support outside the ground space of H top (cf. Fig. 9b). Figure 10: The overlaps | Ψ ± a,b (T )|ψ R | 2 between the final states Ψ ± a,b (T ) of Hamiltonian interpolation and the reference state ψ R (cf. (52)). Observe that the same reference state is used in both figures even though ψ R is naturally associated with the centerpoint in Fig. 10b. The total evolution time is T = 120 in both cases. Comparing with Figs. 9a and 9b, we conclude that throughout the region where the ground space of H top is reached, approximately same state is prepared.
space of H top ): indeed, the latter operator commutes with both H − triv (0, 0) and H top . We also point out that, similarly to the toric code, the local Z i -operators correspond to a combination of pair creation, hopping and pair annihilation of (s, s) anyons.
The preparation stability of the reference state ψ R with respect to the initial Hamiltonians H ± triv (a, b) with negative and positive Z field component is illustrated in Fig. 10. For negative Z field (Fig. 10b) the resulting state Ψ − a,b (T ) has large overlap with the reference state ψ R for almost the entire parameter range. Even when starting from initial Hamiltonians with positive Z field component (Fig. 10a), where the final state does not have a large overlap with the topological ground space (see Fig. 9a), the ground space contribution comes almost exclusively from the reference state. Thus, for doubled semion model, we identify a single stable final state ψ R corresponding to the initial Hamiltonian H triv = − i Z i .

The doubled Fibonacci model
As our last case study of Hamiltonian interpolation, we consider the doubled Fibonacci model described in Section 6.3.4.
(Non)-adiabaticity. Fig. 11 shows the adiabaticity error adia (T ) as a function of the total evolution time T for the initial Hamiltonians H ± triv = ± j Z j . Note that to achieve the same error, the total evolution time T needs to be much longer compared to the toric code and the doubled semion models. It also illustrates that an error of around adia (T ) ≈ 10 −3 is obtained for T = 320: the final state Ψ(T ) overlaps well with the ground space of H top .
In Fig. 12, we consider the non-adiabaticity t → adia (t) along the evolution, again for the initial Hamiltonians H ± triv = ± j Z j . In particular, Fig. 12a, which is for a total evolution time of T = 320, we see that the deviation of the state Ψ(t) from the instantaneous ground state of H(t) can be much larger (compared to the non-adiabaticity adia (T )) along the evolution, even when approaching the end of Hamiltonian interpolation: we have adia (t) 10 −2 for t ≈ 280. The fact that the ground space of the final Hamiltonian H top is reached nevertheless at time t ≈ T is essentially due to the exact degeneracy in the final Hamiltonian H top : In fact, the system is in a state which has a large overlap with the subspace of 'low energy' (corresponding to the 4-fold degenerate subspace of H top ) along the trajectory, but not necessarily with the unique instantaneous ground state of H(t) for t < T . For t = T , the state has a large overlap with the ground space of H top since the latter is higher-dimensional.
This illustrates that the adiabaticity error t → adia (t) along the evolution (i.e., for t < T ) does not provide sufficient information to conclude that the ground space of H top is reached at the end of the evolution. Due to the small energy splitting within the topological "phase" it is more fruitful to view the part of the interpolation close to t ≈ T in terms of degenerate adiabatic perturbation theory [RO14] instead of the traditional adiabatic theorem. Fig. 12a also shows that for T = 320, changing the Trotter time steps ∆t (cf. (48)) from ∆t = 0.1 to ∆t = 0.01 does not significantly change the behavior, particularly for the initial Hamiltonian i Z i . On the other hand, by increasing the Hamiltonian interpolation time for H triv = i Z i to T = 1280, as in Fig. 12b, we see the evolution closely follows the instantaneous ground state. The discrepancy can be seen as a "lag" or delay of the evolved state and the instantaneous ground state and is largest at the "phase transition", H(t) ≈ 1/4H top +3/4 i Z i , where the gap closes.
Note that for this choice of initial Hamiltonians, the vertex terms A v are conserved quantities (as for example in the toric code). Since both |1 ⊗3 and |τ ⊗3 are in the ground space of A v , both signs of the pure Z field lead to a Hamiltonian interpolation which invariantly remains in the ground space of A v . In other words, the adiabaticity error stems from the plaquette terms. Other fields are computationally more costly, since they lift the block decomposition of the interpolating Hamiltonians H(t) induced by the conserved vertex terms, reducing the sparsity of the unitary evolution.
Logical state. Fig. 13 provides information about the final state Ψ ± a,b (T ) of Hamiltonian interpolation, for the family of initial Hamiltonians H ± triv (a, b) (cf. (50)). Again, the figure gives the overlap with a single reference state ψ R . Similarly as before, we choose the latter as the final state of Hamiltonian interpolation, starting with initial Hamiltonian H − triv (0, 0) = − i Z i , and subsequently projected into the ground space and normalized (cf. (52)).
We observe significant overlap of the final state with the reference state ψ R for the whole parameter range for the initial Hamiltonians H − triv (a, b) (Fig. 13b). In contrast, for the initial Hamiltonians H + triv (a, b), the final state depends strongly on the choice of parameters (a, b) (Fig. 13a).
To relate this to the discussion in Section 6 (respectively Conjecture 1), let us first consider the centerpoint of Fig. 13b associated with the initial Hamiltonian H − triv (0, 0) = − j Z j . These terms correspond to a combination of local pair creation, hopping and pair annihilation of (τ , τ ) anyons, as explained in Appendix B. The effective Hamiltonian can be computed at this point based on expression (35) and the S-and T -matrices given in Eq. (45). The result is given numerically in Eq. (100) in the appendix. Computing the ground state ψ eff of this effective Hamiltonian, we observe that with respect to the projections {P 1,1 , P τ,τ , P 1,τ , P τ,1 }, the expectation values of the reference state ψ R and ψ eff are similar, For the whole range of parameters (a, b), the adiabaticity error is small, adia (T ) ≤ 0.005 apart from the point on the boundary of the plot.
The Hamiltonian interpolation reaches the reference state ψ R essentially for the full parameter range.  H top ± 0.001H ± triv (a, b), as illustrated in Fig. 14 (the latter providing an approximate notion of effective Hamiltonians). The figure illustrates that these perturbed states have, as expected, a certain degree of stability with respect to the parameters (a, b). Comparison with Fig. 13 thus points to a certain discrepancy between the behavior of perturbed states and states obtained by Hamiltonian interpolation: Fig. 13a shows high sensitivity of the final state to initial parameters (a, b) (which is absent in the perturbative prediction), whereas Fig. 13b shows that the final state is close to the reference state ψ R throughout (as opposed to the perturbative prediction, where this is not the case along the boundary). To rule out that this discrepancy stems from an insufficiently large choice of the total evolution time T , we also show that different choices of the total evolution time T do not significantly affect the overlap with the reference state along the line b = 0, see Fig. 15.
In summary, we conclude that while for a large parameter range of initial parameters the reference state ψ R is indeed reached, the stability property is less pronounced than for the toric code and the doubled semion models. In addition, a naïve comparison with ground states of perturbed Hamiltonians suggests that the description via effective Hamiltonians does not capture all relevant features. We conjecture that higher orders in perturbation theory are needed to provide more information in the case of the Fibonacci model: the state may be "locked" in eigenstates of such higher-order Hamiltonians before the lowest order effective Hamiltonian dominates.  Fig. 14a, for different values of the total evolution time T . We only compute T = 960 on the rightmost region to show that increasing the total evolution time do not significantly change the final states. It also gives the corresponding overlap | Ψ pert (θ)|ψ R | 2 between the ground state of H pert (θ) and the reference state. The figure illustrates that increasing the evolution time T does not significantly alter the overlap with the reference state. A Equivalence of the self-energy-and Schrieffer-Wolff methods for topological order As discussed in Section 3.4, here we show that at lowest non-trivial order, the expressions obtained from the self-energy-method and the Schrieffer-Wolff method coincide if the Hamiltonian and perturbation satisfies a certain topological order condition. We begin with a review of the exact Schrieffer-Wolff transformation (Section A.1), as well as the expressions resulting from the Schrieffer-Wolff perturbative expansion (Section A.2). In Section A.3, we present some preliminary computations. In Section A.4, we introduce the topological order constraint and establish our main result.

A.1 Exact-Schrieffer-Wolff transformation
As mentioned in Section 3, the Schrieffer-Wolff method provides a unitary U such that preserves the ground space P 0 H of H 0 , and can be considered as an effective Hamiltonian. The definition of the unitary is as follows: let P be the projection onto the ground space of the perturbed Hamiltonian H 0 + εV . Defining the reflections the (exact) Schrieffer-Wolff transformation is defined by the "direct rotation" where the square root is defined with a branch cut along the negative real axis. The effective Hamiltonian is then given by A variational characterization (see [BDL11]) of the unitary U (instead of (54)) is often more useful (e.g., for computing the effective Hamiltonian in the case of a two-dimensional ground space, such as for the Majorana chain): we have U = arg min I − U 2 U unitary and U P U † = P 0 , where A 2 = tr(A † A) is the Frobenius norm.

A.2 The perturbative SW expansion
Since the transforming unitary (54), as well as expression (53), are difficult to compute in general, a standard approach is to derive systematic series in the parameter (the perturbation strength). In this section, we summarize the expressions for this explicit perturbative expansion of the Schrieffer-Wolff effective Hamiltonian obtained in [BDL11]. The perturbation is split into diagonal and off-diagonal parts according to where P 0 is the projection onto the ground space of H 0 , and Q 0 = I −P 0 the projection onto the orthogonal complement. Assuming that {|i } i is the eigenbasis of H 0 with energies H|i = E i |i , one introduces the superoperator Then the operators S j are defined recursively as whereŜ and where Ad S (X) = [S, X]. The constants are a m = 2 m βm m! , where β m is the m-th Bernoulli number. Observe thatŜ The q-th order term in the expansion (6) is where b 2n−1 = 2(2 n −1)β 2n (2n)! . Since our main goal is to apply the perturbation theory to topologically ordered (spin) systems, we can try to utilize their properties. In particular, one defining property of such systems is that, if an operator is supported on a topological trivial region, then it acts trivially inside the ground space. A common non-trivial operation in the ground space corresponds to the virtual process of tunneling an anyon around the torus. This property will allow us to simplify the computation when we want to compute the lowest order effective Hamiltonian. In the following subsections, we will show that although S n is defined recursively based on S 1 , . . . , S n−1 , only the first term −L(Ad V d (S n−1 )) on the rhs of (58) would contribute to the lowest order effective Hamiltonian. The intuition behind this claim is that the other term j≥1 a 2j L(Ŝ 2j (V od ) n−1 ) corresponds to virtual processes which go through the ground space → excited space → ground space cycle multiple times (larger than one). It is intuitive that such virtual processes would not happen when we want to consider the lowest order perturbation.

A.3 Some preparatory definitions and properties
Let G(z) = (zI − H 0 ) −1 be the resolvent of the unperturbed Hamiltonian H 0 . Let E 0 be the ground space energy of H 0 . We set i.e., the inverse is taken on the image of Q 0 . Then L can be written as To organize the terms appearing in the perturbative Schrieffer-Wolff expansion, it will be convenient to introduce the following subspaces of operators.
Definition A.1. For each n ∈ N, let Γ(n) be the linear span of operators of the form where for each j = 0, . . . , n, the operator Z j is either one of the projections P 0 or Q 0 , or a positive power of G, i.e., Z j ∈ {P 0 , Q 0 } ∪ {G m |m ∈ N}. Let Γ (n) ⊂ Γ(n) the span of operators of the form (62) which additionally satisfy the condition i.e., Z 0 and Z n are orthogonal.
For later reference, we remark that operators in Γ (n) a linear combinations of certain terms which are off-diagonal with respect the ground space of H 0 (and its orthogonal complement). In particular, any product of an even number of these operators is diagonal.
The first observation is that the summands the effective Hamiltonian (60) have this particular form.
Lemma A.2. We have and S n ∈ Γ (n) for every n ∈ N (64) for all k, m.
Furthermore, inspecting the Definitions (57) and (56), we immediately verify that Similarly, (63) follows directly from the definitions. We first argue that Inserting the definition of V od and L (that is, (61)), we have where we used the fact GQ 0 = Q 0 G = G and that Q 0 , P 0 are orthogonal projections. This proves the claim (64) for n = 1 and, in particular, shows that S 1 ∈ Γ(1). Similarly, for n = 2, using the definition of V d , a straightforward calculation (using (69)) gives (where h.c. denotes the Hermitian conjugate of the previous expression) and thus with (61) We conclude that (64) holds n = 2 and, in particular, S 2 ∈ Γ(2), as claimed (Eq. (68)).
With (67) and (68), we can use the composition law (66) to show inductively that S n ∈ Γ(n) for all n ∈ N .

A.4 Topological-order constraint
In the following, we will assume that P 0 Γ(n)P 0 ⊂ CP 0 for all n < L . (80) which amounts to saying that (H 0 , V ) satisfies the topological order condition with parameter L (see Definition 3.1). In Section A.4.1, we argue that this implies that the effective Hamiltonian is trivial (i.e., proportional to P 0 ) for all orders n < L. In Section A.4.2, we then compute the non-trivial contribution of lowest order.

A.4.1 Triviality of effective Hamiltonian at orders n < L
A simple consequence of Definition A.1 then is the following.
Proof. It is easy to check that because of property (66), the expression T n 1 · · · T n 2k P 0 is contained in P 0 Γ(n)P 0 , where n = 2k j=1 n j . The claim follows immediately. The argument for P 0 T n 1 · · · T n 2k is identical.
Lemma (A.4) suffices to show that the n-th order effective Hamiltonian H (n) eff is trivial (i.e., proportional to P 0 ) for any order n < L (see Theorem 3.2 below).
A.4.2 Computation of the first non-trivial contribution Lemma (A.4) also implies that the first (potentially) non-trivial term is of order L, and given by P 0Ŝ 1 (V od ) L−1 P 0 . Computing this term requires some effort. Let us define the superoperator V d = −L • Ad V d , that is, For later reference, we note that this operator satisfies as an immediate consequence of (71). We also define the operators Then we can rewrite the recursive definition (58) of the operators S n as where we also introduced Similarly to Lemma A.4, we can show the following: Lemma A.5. Suppose that P 0 Γ(n)P 0 ⊂ CP 0 for all n < L. Then for any for all , m satisfying + m − 1 < L.
Proof. By definition (76), B is a linear combination of terms of the form L(Ŝ 2j (V od ) −1 ) with j ≥ 1, which in turn (cf. (59)) is a linear combination of expressions of the form It hence suffices to show that (The proof of the second statement in (78) is identical and omitted here.) By definition of L, the claim is true if Z 0 = P 0 , since in this case the lhs. vanishes as P 0 V od P 0 = 0. Furthermore, for general Z 0 ∈ {Q 0 } ∪ {G k | k ∈ N}, the claim (79) follows if we can show that P 0 (Ad Sn 1 · · · Ad Sn 2j (V od ))Y V od P 0 ∈ CP 0 , i.e., we can omit L from these considerations. This follows by inserting the expression (61) for L.
Lemma A.6. Suppose that P 0 Γ(n)P 0 ⊂ CP 0 for all n < L. Then for any , m satisfying + m − 1 < L, we have In particular, for every q < L, we have for all k = 0, . . . , q − 2. Furthermore, Proof. By definition of V d , the expression V •m d (B ) is a linear combination of terms of the form and each Z j ∈ {P 0 , Q 0 } ∪ {G m |m ∈ N}. Since A L only involves diagonal operators and the number of factors V d is equal to r < L, we have P 0 A L = P 0 A L P 0 ∈ P 0 Γ(r)P 0 ∈ CP 0 . In particular, where we applied Lemma A.5 with Y = A R (note that A R involves m − r factors V d , and + (m − r) − 1 < L by assumption). We conclude that and since P 0 V •m d (B )V od P 0 is a linear combination of such terms, the first identity in (82) follows. The second identity is shown in an analogous manner.
The claim (83) follows by setting m = k + 1 and = q − k, and observing that + m − 1 = q < L.
Finally, consider the claim (84). We have Lemma A.8. Suppose that P 0 Γ(n)P 0 ⊂ CP 0 for all n < L. Then Proof. By definition (59) and the linearity of Ad · , we havê But by definition of A n−1 , we have if n − 2 < L where we again used the linearity of the involved operations in the second step and Lemma A.6 and Lemma A.7 in the last step (with k = 0 and q = n − 2). Similarly, we have (again by Lemma A.6) if n − 2 < L.
On the other hand, for n = L, we have according to Lemma A.4 and Lemma A.9, hence (89) follows.
B On a class of single-qudit operators in the Levin-Wen model In this appendix, we consider the action of certain single-qudit operators and discuss how they affect states in the Levin-Wen model. For simplicity, we will restrict our attention to models where each particle satisfiesā = a, i.e., is its own antiparticle. Similar local operators were previously considered (for example, in [BSS11]). We introduce the operators in Section B.1 and compute the associated effective Hamiltonians in Section B.2

B.1 Definition and algebraic properties of certain local operators
Recall that for each qudit in the Levin-Wen model, there is an orthonormal basis {|a } a∈F indexed by particle labels. For each particle a ∈ F, we define an operator acting diagonally in the orthonormal basis as As an example, consider the Pauli-Z operator defined in Section 6.3.3 for the doubled semion model. Because the S-matrix of the semion model is given by (see e.g., [ Therefore, the Pauli-Z-operator in the doubled Fibonacci model takes the form where I is the identity matrix. We will write O on ground states of the Levin-Wen model, we used the "fattened honeycomb" description of (superpositions) of string-nets: this gives a compact representation of the action of certain operators (see the appendix of [LW05]), as well as a representation of ground states (see [KKR10]). In this picture, states of the many-spin system are expressed as superpositions of string-nets (ribbon-graphs) embedded in a surface where each plaquette is punctured. Coefficients in the computational basis of the qudits can be obtained by a process of "reduction to the lattice", i.e., the application of F -moves, removal of bubbles etc. similar to the discussion in Section 5. Importantly, the order of reduction does not play a role in obtaining these coefficients as a result of MacLane's theorem (see the appendix of [Kit06]). Note, however, that this diagrammatic formalism only makes sense in the subspace H valid ={|ψ | A v |ψ = |ψ for all vertices v} spanned by valid string-net configurations, since otherwise reduction is not well-defined.
This provides a significant simplification for certain computations. For example, application of a plaquette operator B p corresponds -in this terminology -to the insertion of a "vacuum loop" times a factor 1/D. The latter is itself a superposition of strings, where each string of particle type j carries a coefficient a |ψ is an eigenstate of H top with energy 2. Furthermore, for any ground state |ψ , and any edges e 1 , . . . , e n which (pairwise) do not belong to the same plaquette, the state O (e 1 ) a · · · O (en) a |ψ is an eigenstate (with energy 2n) of H top . The case where the edges belong to the same plaquette will be discussed below in Lemma B.2.
Proof. For concreteness, consider the plaquette operator B p "on the left" of the edge (the argument for the other operator is identical). Because |ψ is a ground state, we have B p |ψ = |ψ . Using the graphical calculus (assuming that the state |ψ is expressed as a string-net embedded in the gray lattice), we obtain Using this fact, we can relate certain products of operators to the string-operator F (a,a) (C) associated with the (doubled) anyon (a, a). That is, assume that the edges {e 1 , . . . , e L } cover a topologically non-trivial loop C on the (dual) lattice (e.g., {10, 12} in the 12-qudit torus of Fig. 4). Then we have

B.2 Effective Hamiltonians for translation-invariant perturbation
According to (91) and (92), a translation-invariant perturbation of the form V = j Z j for the doubled semion or Fibonacci models (as considered in Section 7) is, up to a global energy shift and a proportionality constant, equivalent to a perturbation of the form where a = 1 and the sum is over all edges e of the lattice (Here a = s in the doubled semion model and a = τ in the Fibonacci model). We show the following: where c 1 and c 2 are constants, and the sum is over all topologically non-trivial loops C of length L.
Proof. According to Theorem 3.2, the L-th order effective Hamiltonian is proportional to P 0 (V G) L−1 V P 0 = e 1 ,...,e L P 0 O (e 1 ) a GO (e 2 ) a G · · · GO (e L ) a P 0 up to an energy shift. By the topological order constraint, the only summands on the rhs. which can have a non-trivial action on the ground space are those associated with edges {e 1 , . . . , e L } constituting a non-trivial loop C on the (dual) lattice. Note that for such a collection of edges, every plaquette p has at most two edges e j , e k ∈ {e 1 , . . . , e L } as its sides, a fact we will use below. Our claim follows if we show that for any such collection of edges, we have P 0 O (e 1 ) a GO (e 2 ) a G · · · GO (e L ) a P 0 = cF (a,a) (C) for some constant c.
We show (97) by showing that the resolvent operators G only contribute a global factor; the claim then follows from (94). The reason is that the local operators O such that P 0 O (e 1 ) a GO (e 2 ) a G · · · GO (e L ) a P 0 = Λ k Γ k for k = 1, . . . , L − 1 .
Let |ψ be a ground state of the Levin-Wen model H 0 . We claim that for every k = 1, . . . , L−1, there is a set of plaquettes P k and a constant c k (independent of the chosen ground state) such that |ψ .
(ii) The (unnormalized) state |ψ is an eigenstate of H 0 . Its energy k is independent of the state |ψ .
(iii) The set P k only contains plaquettes which have two edges in common with {e k , . . . , e L }.
Note that for k = 1, this implies P 0 O

71
Assume now that (i), (ii) hold for some k ∈ {2, L}. Then we have according to (98) where It hence suffices to show that for some choice of plaquettes P k−1 , we have |ψ is an eigenstate of H 0 with energy k−1 (independent of |ψ ).
(c) that the set P k−1 only contains plaquettes sharing two edges with {e k−1 , . . . , e L }.
By assumption (iii) and the particular choice of {e 1 , . . . , e L }, none of the plaquettes p ∈ P k contains the edge e k−1 . Therefore, we can commute the operator O We then consider two cases: • If e k−1 does not lie on the same plaquette as any of the edges {e k , . . . , e L }, then application of O (e k−1 ) a creates a pair of excitations according to Lemma B.1 and the state (99) is an eigenstate of H 0 with energy k−1 = k + 2 > E 0 . In particular, setting P k−1 = P k , properties (a)-(c) follow.
• If there is an edge e , ≥ k such that e k−1 and e belong to the same plaquettep, then the state (99) is a superposition of states with Bp excited/not excited, that is, we have |ϕ = (e k−2 ) a since these do not share an edge withp, hence Λ k−1 (I − Bp)|ϕ = 0 (recall that Λ k−1 = P 0 Λ k−1 includes a projection onto the ground space). Thus setting P k−1 = P k ∪{p}, we can verify that (a)-(c) indeed are satisfied. (The case where there are two such plaquettesp can be treated analogously.) Let us compute the effective Hamiltonian (96) for the case of the rhombic torus, or more specifically, the lattice we use in the numerical simulation, Fig. 4. It has three inequivalent weight-2 loops: {10, 12}, {1, 2}, {5, 7}. Follow the recipe in Section 6.2, respectively Section 7.2, these three loops are related by a 120 • rotation. The corresponding unitary transformation for this rotation is given by the product of matrices A = T S when expressed in the flux basis discussed in Section 6.2 (for the doubled Fibonacci model, the latter two matrices are given by (45)). Similarly, we can express the action of F (a,a) (C) in this basis using (36), getting a matrix F . By (39), the effective Hamiltonian for the perturbation − j Z j is then proportional to (when expressed in the same basis) Note that the overall sign of the effective Hamiltonian is not specified in (34), but can be determined to be negative here by explicit calculation. For example, substituting in the S matrix (Eq. (45)) of the doubled Fibonacci model, we have F = diag(ϕ + 1, −1, −1, ϕ + 1) for the Fibonacci model. It is then straightforward to obtain the ground state of H eff , which is 0.715|(1, 1) + (0.019 − 0.057i)|(τ, 1) + (0.019 + 0.057i)|(1, τ ) + 0.693|(τ, τ ) , (100) where |(a, b) is a flux basis vector, i.e., the image of P (a,b) (C) (see Section 6.2) up to some phase.