Incorporating non-local anyonic statistics into a graph decomposition

In this work we describe how to systematically implement a full graph decomposition to set up a linked-cluster expansion for the topological phase of Kitaev’s toric code in a field. This demands to include the non-local effects mediated by the mutual anyonic statistics of elementary charge and flux excitations. Technically, we describe how to consistently integrate such non-local effects into a hypergraph decomposition for single excitations. The approach is demonstrated for the ground-state energy and the elementary excitation energies of charges and fluxes in the perturbed topological phase


Introduction
The research on topological quantum phases [1][2][3] is a very dynamic field in modern physics.Such exotic phases have highly entangled ground states that cannot be described by Landau spontaneous symmetry breaking and possess anyonic elementary excitations with non-trivial particle statistics [4,5] deviating from conventional bosons and fermions.These anyons are the essential ingredient for the field of topological quantum computing exploiting non-Abelian braiding properties [6,7].Topological order plays a crucial role in the physics of the fractional quantum Hall effect [8,9] and in frustrated quantum systems as quantum spin liquids in quantum simulators or condensed matter realizations [10,11].In this context, strongly correlated Mott insulators with strong spin-orbit interaction like the iridates [12,13] as well as α-RuCl 3 [14][15][16] have been intensively studied in recent years.These quantum materials are discussed as potential realizations of Kitaev's honeycomb model [17], which is exactly solvable possessing topologically ordered ground states.In the limit of strong anisotropic interactions Kitaev's honeycomb model can be mapped to the well-known toric code [6] arising as an effective low-energy model in perturbation theory.
One prominent method which has been applied successfully to determine the quantumcritical properties of the 2D toric code in the presence of a uniform magnetic field [21,22,25] are high-order series expansions applying the method of perturbative continuous unitary transformations [47,48].In these works the excitation energies of elementary charges and fluxes have been calculated as a series in the low-field limit.The gap closing is then analyzed to extract quantum critical points and associated critical exponents.The non-local effects of the braiding statistics have been taken into account in a post-processing procedure by correctly incorporating the winding of charges around fluxes (or vice versa).This has been done with large clusters, e.g., using Entings finite-lattice method [49], in order to define particle strings in a fixed gauge-invariant way so that the particle as well as the relevant fluctuations are both present on each cluster.
However, linked-cluster expansions are often implemented most efficiently by applying a full graph-decomposition so that calculations on individual graphs cost minimal memory and time.In this work we describe how to integrate systematically the non-local effects of the anyonic statistics in a hypergraph decomposition.The latter have recently been established to naturally allow the treatment of multi-site perturbations within linked-cluster expansions [50].Interestingly, they are also suited to set up a linked-cluster expansion to determine the energy of single charge and flux excitations for the topological phase of Kitaev's toric code in the presence of a general uniform field.Specifically, we determine the ground-state energy and the elementary excitation energies of charges and fluxes in the perturbed topological phase.
The paper is organized as follows.In Sec. 2 we introduce the toric code in a field including an extensive introduction to the bare toric code.In Sec. 3 we present all generic aspects for series expansions of the perturbed topological phase in the toric code at finite magnetic fields while Sec. 4 contains all information for a linked-cluster expansion with full hypergraph decomposition.Final conclusions are drawn in Sec. 5.

Kitaev's toric code in a field
We investigate the perturbed topological phase of the toric code in a uniform field.The Hamiltonian of this model is therefore the sum of the toric code H TC and a general uniform magnetic field where ⃗ h = (h x , h y , h z ) T ∈ R 3 and ⃗ σ = (σ x , σ y , σ z ) T .In contrast to the bare toric code described in Subsec.2.1, this model is not exactly solvable and displays a rich quantum phase diagram [25].For the special case of a magnetic field pointing in xor z-direction the model is isospectral to the well-known transverse field Ising model on the dual square lattice in the relevant low-energy sector [18,19] displaying a second-order quantum phase transition in the 3D Ising universality class.The quantum phase diagram of the toric code in a magnetic field in the xz-plane is obtained with series expansion methods [21], quantum Monte Carlo simulations [23,24], and tensor networks [26].Notably, the universality class of the quantum phase transition remains 3D Ising except for the symmetric case h x = h z where a multicritical point is expected [21,23,25].In contrast, for a magnetic field pointing in transverse y-direction, the model is dual to the Xu-Moore model [51] as well as to the quantum compass model [52] and the nature of the phase transition is strongly first order [22].The full extent of the topological phase in the presence of a general uniform magnetic field has been determined by combining high-order series expansions and tensor network algorithms [25] displaying rich physical behavior depending on the magnetic field direction.This includes planes of first-and second-order quantum phase transitions as well as multicritical lines.
In this work we focus on the perturbed topological phase at finite fields and we describe how to set up high-order series expansions.Consequently, we describe the non-trivial unperturbed limit ⃗ h = 0 of the bare toric code next.

Kitaev's toric code
Kitaev's toric code is defined on a square lattice with spin-1/2 degrees of freedom on the edges of the lattice [6].The Hamiltonian is given by The star operators X and the plaquette operators Z are both defined as products of four Pauli matrices where σ x i , σ z i are the usual Pauli matrices.The first product runs over the four sites surrounding a vertex, the second one over the four sites at the edges of a plaquette as illustrated in Fig. 1.As these operators are products of Pauli matrices their eigenvalues x , z are equal to ±1 [6].Furthermore, all these operators commute mutually As a consequence, one can construct a state where all eigenvalues x , z are equal to +1, which corresponds to a ground state where N is a normalization factor and |ref〉 is a reference state.This reference state has to be chosen such that it is not orthogonal to the ground-state space.For example, the states |⇑〉 or |⇒〉, where all spins point in positive zor x-direction can be used.The ground state |GS〉 is unique on an open plane.However, the ground-state manifold features a non-trivial topological degeneracy which equals 4 g on a compact orientable surface of genus g indicating topological order [6].Importantly, products of operators X and Z act trivially on the ground state.Furthermore, any contractible loop of σ x or σ z matrices is equal to the product of operators X or Z contained in the loop respectively [6].The ground-state energy is given by E tc 0 = −N /2 with N the number of spins which equals the total number of plaquettes and stars.
Acting with σ x i (σ z i ) on site i of the ground state creates two flux (charge) excitations corresponding to z = −1 (x = −1) adjacent to site i.These topological excitations are located at the plaquettes (vertices) of the lattice and behave like hardcore bosons with a mutual anyonic statistics [6,25,35].More generally, acting with string operators on the ground state, where p (p) are open paths on the (dual) lattice, creates two excitations at the end of the respective path [6].Examples for two such operators are illustrated in the right panel of Fig. 1.Note that it is not possible to have an odd number of excitations on the torus [6].However, on an open plane, one can create a pair of excitations and move one of them to infinity.A single charge x = −1 at position ⃗ r is denoted by the state |⃗ r, 〉.This state is uniquely defined up to the concrete operator which creates the excitation from the ground state.In order to define such an excited state unambiguously we define it by a sequence of Pauli-operators acting on the reference ground state where p ⃗ r is a straight open path which goes from the excitation to infinity in negative xdirection as illustrated in the left panel of Fig. 2.This way we fixed the gauge freedom to define these states uniquely.In the following we call these states canonical one-charge states.
In the same fashion one can define canonical one-flux states where p⃗ r is a straight open path on the dual lattice which is going straight into negative xdirection from the position ⃗ r of the flux excitation to infinity as illustrated in the right panel of Fig. 2. We stress that non-canonical one-particle states with a different string can always be For a non-canonical one-charge state which is obtained by acting with σ z i along an open path l ⃗ r ̸ = p ⃗ r from infinity to ⃗ r this can be achieved in the following way where Γ z is chosen such that SciPost Physics

Submission
One specific example is illustrated for the case of a single-charge state in Fig. 3.The analogue procedure can be done for corresponding single-flux states.More generally, applying the same procedure and conventions, one can generate multiparticle states with an arbitrary number of charges and fluxes.This can be done by successively acting with the corresponding (canonical) string operators on the ground state |GS〉.
Next we introduce a specific representation of these canonical multi-particle basis states which is convenient for high-order series expansions discussed in the following sections.To this end we consider as reference states |ref〉 the fully polarized states |⇑〉 and |⇒〉 in zand x-direction, respectively.Once the reference state is chosen, all canonical states are uniquely defined including the overall phase.Let us consider an arbitrary multi-particle state with energy only depending on the eigenvalues of plaquettes and stars.In our construction these states can be written as the product of (1± X ) (1± Z ) and the action of string operators (for each particle one) on the reference state |ref〉.The second part of this product corresponds to a spin product state which we call spin background.It can be described by the set of spin eigenvalues {s i } with s i ∈ {±1} on all sites i.For the specific case of a two-particle state with one charge and one flux, one has where 1 is the star at position ⃗ r 1 , 2 is the plaquette located at ⃗ r 2 , and σ x i |ref〉 is the spin background.In general, we keep track of the three sets {x }, {z }, and {s i } of 2 degrees of freedom.This way we over-parametrize the states.However, the spin background will turn out to be useful in the next sections to treat the non-trivial statistics of charges and fluxes within the high-order series expansions.

Series expansions in the perturbed topological phase
Our goal is to describe a systematic way how to incorporate non-local anyonic statistics into a high-order series expansion exploiting a full graph decomposition considering the topological phase of the toric code as the unperturbed starting point.How high-order series expansions for arbitrary field directions can be derived using finite clusters is explained in [53], where also Entings finite lattice method [49,54] has been applied.Such series expansions have been applied successfully to the toric code in a field [21,22,25].However, to the best of our knowledge no full graph decomposition has been performed for general field directions.In this section we review how to calculate the ground-state energy and one-quasiparticle excitation energies in the topological phase using perturbation theory about the low-field limit.To this end we consider perturbative calculations on finite clusters which are sufficiently large to contain the relevant physical properties in the thermodynamic limit [48,54,55].The discussion of the full graph decomposition of the perturbed topological phase is then presented in Sec. 4.

Set up
We perform perturbative calculations about the limit ⃗ h = 0. We therefore take the canonical states |{x }, {z }, {s i }〉 with energy E tc ({x }, {z }) as unperturbed basis (see Sec. 2.1 for details).We stress again that the unperturbed energies of these basis states do not depend on the spin background {s i }.The perturbation V ≡ − i ⃗ h • ⃗ σ i is a sum over local terms acting on individual sites i. Locally, besides the obvious action on the spin at site i, the Pauli matrix σ x i (σ z i ) flips the two eigenvalues z (x ) attached to site i while the Pauli matrix σ y i flips all four eigenvalues of star and plaquette operators containing site i.We can therefore split the perturbation V as follows so that [H TC , T n ] = nT n and n corresponding to the net change of the total charge and flux particle number due to the action of T n .Consequently, at finite fields, the system becomes a challenging quantum many-body problem where the bare number of charges and fluxes is not conserved and the elementary charge and flux excitations gain a finite dispersion and interact.The properties of the unperturbed toric code as well as the decomposition Eq. 14 meet all the requirements to apply the method of perturbative continuous transformation (pCUTs) [47,48], which has been done in a series of works [21,22,25].The pCUT method allows to map (1) to an effective quasi-particle-number conserving Hamiltonian H eff perturbatively exact up to the calculated order in the parameters h x , h y , and h z .This effective Hamiltonian obeys [H TC , H eff ] = 0, i.e., the effective Hamiltonian is block-diagonal and the quantum manybody problem reduces to a few-body problem in terms of dressed charge and flux excitations.Specifically, H eff can be written as follows where m ≡ (m 1 , . . ., m k ) with m i ∈ {0, ±2, ±4}, T (m) = T m 1 T m 2 T m 3 . . .T m k is a product of k operators T n in order k perturbation theory, M (m) = m i , and C(m) rational numbers.In each perturbative order k one therefore has a weighted sum of quantum fluctuations described by the perturbative process T (m).
The Hamiltonian (15) is not normal ordered so that physical properties can not be extracted directly.The normal ordering is most efficiently done by calculating matrix elements on finite clusters as only linked processes can have a non-vanishing contribution to these matrix elements.Indeed, if a finite cluster is sufficiently large so that all linked processes of a given perturbative order fit on this cluster for a specific matrix element, then this matrix element can be calculated and is directly valid in the thermodynamic limit.This therefore allows to obtain the normal-ordered form of the effective Hamiltonian.
The appearance of linked processes can be directly seen by introducing where τ n,b (α) acts on bonds b (α) with α ∈ {x, y, z}.These different bond types are illustrated in Fig. 4. The bond b (α) is the ordered set of charge sites , flux sites , and spin sites i on which the operator τ n,b (α) acts.Clearly, any bond contains the spin site i on which the respective Pauli matrix acts.Additionally, all charge and flux sites hosting operators whose eigenvalues x , z are affected by the respective Pauli matrix need to be included in the bond.Note that it is important to distinguish the different types of sites within the bonds.Here it is necessary to distinguish spin sites from flux and charge sites, as the operators act differently on them.Accordingly, a bond b (z) (b (x) ) includes a spin site and the neighboring charge sites (flux sites), whereas a bond b ( y) contains a spin site and the four neighboring charge and flux sites.The effective Hamiltonian can then be exactly rewritten as In order to determine prefactors of the normal-ordered quasi-particle conserving operators in H eff via the calculation of matrix elements, one has to evaluate where |Φ〉 and |Ψ〉 are states on sufficiently large clusters with the same number of quasiparticles.One therefore has to specify the action of the local operators τ n,b (α) on the canonical states |{x }, {z }, {s i }〉.We stress that the action on the spin background depends on the choice of the specific reference state |ref〉.As an example, we list in Tab. 1 the different processes taking the fully polarized state |⇒〉 as the reference state.
As outlined in Sec. 2, for the perturbed topological phase of the toric code in a field, one has to take care that scalar products are made with canonical states.States arising from the action of the operators τ n,b (α) are not necessarily canonical.In the following subsections we discuss how to evaluate these scalar products properly for the ground-state energy and for one-quasi-particle states on a single large cluster.In principle the scalar products for higher quasi-particle states are based on the same principles and are described in the literature [53].

Ground-state energy
Let us start the discussion for the ground-state energy E 0 by focusing on a single term of the effective Hamiltonian H eff in order k perturbation theory.Ignoring prefactors in Eq.( 15), this demands to calculate the following expectation value where |GS〉 represents the unique ground state of the toric code on the chosen open cluster.
Table 1: On the left the action of σ x i on the local state configuration |z 1 , z 2 ; si 〉 on a b (x) -bond, which contains two plaquette operators and the background spin s i .On the right the action of σ z on the local state |x 1 , x 2 ; s i 〉 on a b (z) -bond type is given, which contains two star operators and the same background spin.The action on the spin background variables s i depends on the reference states |ref〉.Here we have used the fully polarized states |⇒〉 in x-direction.The upper, middle, and lower line refers to the processes of τ 2,b (α) , τ 0,b (α) , and τ −2,b (α) with α ∈ {x, z}, respectively.For the sake of clarity we use Boolean values z = 1 2 (1 − z ) and x = 1 2 (1 − x ), which count the number of excitations at a given or and for the background spins si = 1 2 (1 + s i ).Furthermore, σ y = iσ x σ z acts on the background spin s i as well as on the four eigenvalues of the star and plaquette operators which surround site i.
The successive action of the τ-operators in Eq. ( 17) is properly treated on the state level by the processes listed in Tab. 1.Because the product of τ-operators is quasi-particle conserving, it can not create any excitations.Instead, the final state |f〉 ≡ τ m 1 ,b (α 1 ) either zero or is the ground state |GS〉 with an additional phase factor [i n a ] with n a ∈ {0, 1, 2, 3} resulting from the action of the τ-operators.
One can explicitly write down the final state |f〉 as The non-contributing sequences are easily identified, so we restrict the discussion to the contributing ones.For the explicit reference state |⇒〉, the action of σ on the reference state is given as |s x 〉 is a product state in the σ x -basis and [i n a ] with n a ∈ {0, 1, 2, 3} is the phase factor resulting from the action of the Pauli matrices, whereas Γ z is a set of plaquettes such that Note that we are in the ground-state sector, so |s x 〉 can at most contain several contractible loops of flipped spins in the spin-background, which can be replaced by a product of operators Z .Inserting Eq. ( 19) into Eq.( 18) yields Using the operators defined in Tab. 1 results exactly in the middle part of Eq. ( 21), i.e. they yield the result in terms of the state ∈Γ z Z |GS〉 and correctly include the additional phase [i n a ] within the amplitude of this state.
Interestingly, the state |f〉 can also be evaluated for an arbitrary reference state.To this end we sort the product of Pauli matrices in Eq. ( 18) by their flavor taking into account the nontrivial commutation relations.Note that Pauli matrices σ y can always be written as iσ x σ z .In the ground-state sector one can always write the obtained sorted expression as the product of contractible loops of σ z or σ x .Each contractible loop operator can be replaced exactly by the product of plaquette or star operators contained in the respective loop Since any product of star and plaquette operators yields identity on any ground state, one directly gets the canonical ground state |GS〉 with an additional phase factor [i n a ] with n a ∈ {0, 1, 2, 3}, which stems from the sorting of the Pauli matrices So once the phase factor is determined, no post-processing is necessary in the ground-state sector and one has 〈GS| τ m 1 ,b (α 1 )

One-quasi-particle energies
The one-quasi-particle sector of the effective Hamiltonian H eff represents two one-particle problems: One for a single charge quasi-particle and one for a single flux quasi-particle.Both problems are exactly decoupled since the magnetic field does not contain any process transforming a charge into a flux or vice versa due to parity conservation of fluxes and charges.Furthermore, there is an exact self-duality in the toric code in a field so that charge and flux excitation energies are identical up to an interchange of h x and h z .Let |⃗ r 1 〉 be either the canonical one-charge state |⃗ r 1 , 〉 or the canonical one-flux state |⃗ r 1 , 〉.In real space one then has to calculate the one-particle hopping amplitudes a ⃗ δ given by a Introducing corresponding one-quasi-particle states in Fourier space with N the number of spins, the one-particle dispersion of charges and fluxes is The minimum of this dispersion is called the one-particle gap ∆ located at momentum ⃗ k = 0 for charges and fluxes.In the following we call the charge gap ∆ and the flux gap ∆ .Next we discuss how to evaluate the matrix elements a ⃗ δ .As in the calculation of the ground-state energy, we can focus on the contribution of a single perturbative term to the hopping amplitude a ⃗ δ with ⃗ r 2 ≡ ⃗ r 1 + ⃗ δ.In order to evaluate the scalar product 〈⃗ r 2 |f〉, one has to express |f〉 in terms of a canonical one-quasi-particle state along the lines discussed in Sect. 2.

First let us rewrite a canonical one-charge state as
For a given reference state |ref〉 this defines a canonical spin-background for one-charge states Next we investigate the action of a non-vanishing quasi particle number conserving operator sequence on a canonical one-charge state where 2 is the charge site at position ⃗ r 2 .While it is clear that the charge moved from ⃗ r 1 to ⃗ r 2 , the spin-background is not necessarily in a canonical form.At this point it is most instructive to use |⇒〉 as reference state and to bring the spinbackground into a canonical form.For this choice we have Using σ y = iσ x σ z , the action of the Pauli matrices on |⇒〉 is evaluated where |s x 〉 is a product state in the σ x -basis and the phase factor i n a with n a ∈ {0, 1, 2, 3} results from the action of the Pauli matrices.
The expression on the right results from the fact that any |s x 〉 can be written as a product of σ z matrices acting on the reference state |⇒〉.In the current case this product of σ z can always be chosen as the canonical Pauli string operator on p ⃗ r 2 of the one-charge state |⃗ r 2 , 〉 times some product of Z over the appropriate set Γ z .With this the final state |f〉 becomes Note that the result is again obtained in terms of a canonical state |⃗ r 2 , 〉 and an additional phase factor [i n a ] stemming from the action of the Pauli matrices on the reference state as detailed in Eq. 30.The operators given in Tab. 1 yield the result in terms of a state | f〉 with a phase So although the spin-background of | f〉 is not necessarily canonical we have reflecting that canonicalizing the state | f〉 is a trivial operation due to the absence of fluxes.Note that this is different if we choose the reference state |⇑〉 and act with the Pauli matrices to the right Then we can write where the factor [i n a ] comes from the action of the Pauli operators.We can then commute the product over X to the left which is a product of X acting on a one-charge state.And hence we get the following expression where Clearly, in this convention, the definition of the operators in Tab. 1 needs to be adapted.These adapted operators then also yield the result in the form However, in contrast to Eq. ( 33), we have So we still need to determine n b from the spin-background of the state | f〉 by checking whether there are contractible loops surrounding the charge at ⃗ r 2 .Physically, each such loop corresponds to the winding of a flux around the charge.As a consequence, the total effect of all these loops depends only on the parity of the loop number.In case of even (odd) parity no (an) additional sign results which has to be included in this non-trivial postprocessing procedure.
In the previous part we used an explicit reference state because it is easy to act with sequences of Pauli operators on product states in the spin basis.But actually the arguments are independent of the reference state.Note that we effectively just sorted the Pauli operators by flavor.Splitting σ y = iσ x σ z we arranged the product σ such that σ x i acts before σ z i or vice versa.Accordingly, on an operator level, the signs stem from the ordering of the Pauli operators.In this sense the τ-operators implement a bookkeeping technique for the necessary reordering and yield the appropriate sign for the ordered product.
As for the ground-state energy, the state |f〉 can then be written independently of the reference state as where ς α is the part consisting of σ α .The additional factor [i n a ] stems from the commutation relations of the Pauli matrices.It is clear that ς x corresponds to contractible loop operators and hence we have The product ς z corresponds to a product of stabilizers over the set Γ z with a canonical path going from ⃗ r 2 to infinity.One therefore has As a next step we move the product of X and Z to the left where n b captures whether ∈Γ x X and i∈p ⃗ r 2 σ z i commute or anticommute What remains is the action of stabilizer operators on a canonical one-charge state where As n b ∈ {0, 1} we have which means that the phase factor stems from the sorting of the Pauli matrices.Note that we did choose a convention on the sorting of the Pauli matrices in this calculation, by sorting σ x to the right.The operators given in Tab. 1 take into account the signs arising from this reordering and yield the final state in the form Comparing to Eq. ( 48), one sees that the canonicalization of | f〉 cannot yield an additional phase factor.Overall, these considerations allow to set up series expansions of the perturbed topological phase in the zero-and one-particle sector using sufficiently large clusters.A generalization to many particles requires further considerations treating the spin background properly in the evaluation of scalar products between multi-particle states [35,53].Next we go beyond series expansions on large clusters by performing a full graph decomposition.

Hypergraph decomposition for the perturbed topological phase
In this section we describe how to execute linked-cluster expansions [56,57] for the perturbed toric code in the topological phase at finite fields using a full hypergraph decomposition [50].In this work we restrict ourselves to the zero-and one-particle sector.In these sectors it turns out that the non-trival mutual statistics can be correctly included by taking the semi-infinite string operators into account within the hypergraph decomposition.

Linked-cluster expansions
On a sufficiently large cluster C as introduced in Sec. 3, the ground-state energy as well as the irreducible one-particle matrix elements of the effective Hamiltonian can be written as [56] where κ(C) is the reduced contribution of the cluster C, i.e., the contribution which contains only the perturbative processes which act at least once on every bond of the cluster C [56,58].
Note that the sum runs over all sub-clusters c of C and not over all distinct sub-clusters of C.
For the ground-state energy, the contribution K(c) of a cluster c is where H c eff ist the restriction of H eff to the bonds and sites of the cluster c.Accordingly, |GS〉 is the unperturbed ground state on the full system.Obviously, H c eff can also be evaluated on the restriction of |GS〉 to the cluster c.For the one-quasiparticle states, we subtract the groundstate energy and calculate matrix elements of the form [59] The contribution K(c) of cluster c is then where |⃗ r〉 is a canonical state with a flux or charge at position ⃗ r defined on the full system.The reduced contribution κ(C) can be calculated by where the sum runs over all proper sub-clusters of c, which means that the reduced contribution depends only on reduced contributions of smaller clusters and can be determined recursively.The restrictions of the Hamiltonian on isomorphic (structurally equivalent) clusters are equivalent [56] and thus equivalent matrix elements of the (effective) Hamiltonian on these clusters agree.So we can rewrite (50) as where E labels an equivalence class of sub-clusters, N (E, C) gives the number of sub-clusters of C which belong to E, and κ(E) is the reduced contribution of a cluster from the class E [56,57].
For extensive quantities N (E, C) is normalized (usually to the number of sites in the cluster C).Importantly, all clusters c and c ′ which are in the equivalence class E fulfill The distinction of the equivalence classes is based on the structure of the clusters and the states on the clusters within the considered matrix element.To this end the clusters are usually associated to graphs, where the vertices represent sites and bonds linking these sites correspond to edges.As structural equivalent clusters are represented by isomorphic graphs we refer to them as isomorphic clusters.Further, the local states on the clusters within the considered matrix element are incorporated using additional vertex colors.The equivalence of clusters then corresponds to (color-preserving) isomorphism of the respective graphs.
However, for the problem at hand the bonds b (x) , b ( y) , b (z) join multiple sites, and not necessarily in a symmetric way, so there is also some orientation within the bonds.A representation in terms of hypergraphs [60] appears natural for multi-site couplings.Interestingly, hypergraphs can be unambiguously represented by bipartite graphs, which are referred to as the König representation [61,62].So instead of the usual graph representation we use an enhanced König representation to represent the clusters by graphs as described in [50].Also in this case information about the local states on the clusters are incorporated into the graph representation using additional vertex colors.We concretely discuss how to construct the graph representation of the clusters including some concrete examples in Sec 4.3.
Furthermore, for a given perturbation order, the sums in Eqs. ( 50) and ( 55) need to consider only clusters up to a given number of bonds, because a cluster can not contribute in perturbation orders smaller than its number of bonds [56].As the pCUT method is clusteradditive, κ(C) vanishes for any disconnected cluster C and so it is sufficient if these sums run over all relevant connected sub-clusters or the corresponding equivalence classes respectively.Typically the cluster C is chosen large enough to host all perturbative contributions in the desired perturbation order, so that the factors N (E, C) on the finite cluster C are the same as on the infinite lattice L (after appropriate normalization in the case of extensive quantities).Additional selection rules [57,63,64] to check whether a cluster (or an entire equivalence class) can contribute at a given order are important to further truncate the sums and increase the efficiency of the method.We give some details about such heuristics for the toric code in a field in App. A.

Contributions from individual clusters
For the low-field limit of the toric code in a general uniform field the unperturbed Hamiltonian depends on the eigenvalues x , z which are defined locally at the charge or flux sites.The general uniform field then acts as a perturbation on three different types of bonds as described in the last section.These three bonds correspond to the elementary building blocks of the clusters.This means every cluster relevant for a perturbative expansion is determined by a set of bonds, which is a subset of all the bonds in the lattice.Clusters consisting of a single site are relevant only in zeroth order, so their contributions can be easily calculated and are neglected in the following discussion.Note that in contrast to the hypergraph decomposition in Ref. [50], we have several types of sites in the current problem: The spin-sites s i , the flux-sites , and the charge-sites .
In Sec. 3 we describe how to calculate the ground-state energy and one-quasiparticle matrix elements on finite clusters.In principle each of these calculations can be decomposed as given in Eq. ( 50).However, for a given irreducible matrix element, the contribution of a cluster can depend on particles residing outside of the clusters in the current case, which is counterintuitive and a consequence of the non-local properties of topological phases.Such effects are however absent in the ground-state sector.Indeed, we showed that for the ground-state energy it is sufficient to act with the operators τ n,b (α) (see Tab. 1) to the right so that the result is given in terms of the canonical ground state |GS〉 times an additional phase factor.So the ground-state energy can easily be evaluated on each individual cluster and does not depend on any degree of freedom which is outside of the respective cluster.
In the one-quasiparticle sector, the procedure depends on the choice of the reference state.First, there are the specific reference states where the strings of the considered quasi-particle are visible but the corresponding strings of the other quasi-particle type are not visible in the spin-background.Let us consider the generic example of a single charge excitation taking the reference state |⇒〉.The phase factor [i n a ] in the one-charge state in Eq. ( 31) originates from the action of the local operators τ n,b (α) on the canonical state.Due to the absence of fluxes and the specific choice of the reference state, this one-charge state is a canonical state up to the phase [i n a ].As a consequence, we can evaluate these irreducible one-charge hopping elements on individual clusters without the necessity to take degrees of freedom outside the cluster into account.
In contrast, using the reference state |⇑〉 as in Eq. ( 36), the phase factor [i n a +2n b ] stems from the evaluation of the τ n,b (α) on this reference state and the canonicalization of the resulting state.In this case the state resulting from the action of the τ n,b (α) may contain loops which surround the charge and hence lead to an additional phase factor of −1.In the formulation for arbitrary reference states, these conventions correspond to different sorting of the product of Pauli matrices acting on the reference state.In the first convention σ x is sorted to the right of σ z , in the second convention it is the other way around.Overall, in the appropriate convention the evaluation of the sequence of τ n,b (α) yields a state which is equivalent to a canonical state and an additional phase factor.So these irreducible hopping elements are easily evaluated on each cluster and do not depend on any degree of freedom outside of the cluster.Clearly, by the duality of charges and fluxes also one-flux irreducible matrix elements can be determined on the individual clusters, by choosing the appropriate convention.
Note, however, that in order to evaluate the contribution of a cluster we need to take into account all the degrees of freedom within the cluster including the spin-background variables.
Especially for the irreducible one-quasiparticle hopping elements this leads to the observation, that clusters which do not host particles contribute due to a non-trivial spin-background on the cluster.Due to this spin-background the respective matrix element of the effective Hamiltonian might evaluate differently to the ground-state energy (which is subtracted to obtain the irreducible element).A physically simple, but enlightening example is a cluster which −1 Figure 5: Example for a perturbation theory process which contributes to the excitation gap, although it does not directly act on the quasiparticle.The process consists of twelve σ x operators which act on the red sites forming a loop around the excitation and interacting with the string operator.This operator consists of σ z operators acting on the blue sites.Both operators, the loop and the string operator, act on the site which is red and blue.winds around a charge as illustrated in Fig. 5.It is clear, that a process winding a flux around the charge results in a different contribution, than the same process without the charge even though the process does not directly involve the respective charge site.So it is consistent that these clusters contribute to the irreducible hopping elements.Note that this type of process already occurs in order four perturbation theory.

Identification of equivalent clusters
In the previous subsection we have shown that we can indeed write an expansion in the form of Eq. ( 50).We also explained that for the ground-state energy and the one-quasiparticle matrix elements we can always choose a reference state, such that we can evaluate the contributions of the individual clusters without taking into account degrees of freedom which are not on the cluster.Throughout this subsection we assume to work in such a convention.The next step is to exploit symmetries and structural equivalence to identify clusters which have the same contribution to the desired matrix element to arrive at an expansion of the form of Eq. ( 55).The distinction of clusters is based on the hypergraph decomposition explained in Ref. [50].In contrast to the hypergraph decomposition in Ref. [50], we have several types of sites in the present problem: spin-sites s i , flux-sites , and charge-sites .From a technical point of view it is interesting that a given perturbation operator can act differently on different types of degrees of freedom, namely the spin-background variables and the eigenvalues of the star and plaquette operators.So while a cluster still consists of sites and bonds, the bonds have some orientation.This means different sites have different roles in the bonds.While the eigenvalues of the stabilizer operators are always flipped, the action on the spin-background can be different.So the structure of the bonds cannot be captured by ordinary hyperedges,  which are just subsets of the vertex set, in a hypergraph representation of the cluster.Interestingly, several definitions of oriented hypergraphs can be found in the literature [65][66][67][68][69][70][71].
Here, we assign labels to incident edge-vertex pairs, to include the roles of the vertices within the bonds.Hypergraphs can be represented by bipartite graphs, which are also referred to as the König representation of a hypergraph [61,62].The two parts of the bipartite König graph represent the edges (edge-part) and the vertices (vertex-part) of the hypergraph.Two vertices in the König representation are adjacent if and only if they correspond to an incident edge-vertex pair within the hypergraph.Importantly, two hypergraphs are isomorphic if and only if their König representations are isomorphic [61,62].
The labels for the incident hyperedge-vertex pairs can be incorporated using edge colors within the König representation.Note that this is consistent with the fact that the edges in the König representation correspond to incidences of edge-vertex pairs in the original hypergraph.Furthermore, we use vertex colors to distinguish the two parts of the bipartite graph and to incorporate the different bond-types into the edge-part of the König representation.
As an example we give the König representation of the clusters which consist of a single bond in Fig. 6.In the end this enhanced König representation is the basis for the distinction of equivalence classes of clusters within the full hypergraph decomposition presented here.For calculations of the ground-state energy this representation of the clusters is indeed sufficient.
For the hopping elements we need to explicitly incorporate the information on the involved states on the clusters, as the contribution of a clusters also depends on the states considered within the matrix element.Finally also this can be incorporated via vertex colors.

Ground-state energy
For the ground-state energy equivalent clusters are identified using hypergraph isomorphism accounting for the different roles of the vertices within the hyperedges and the different bond types b (α) with α ∈ {x, y, z}.In terms of the König representation this means that the isomorphism preserves the edge and the vertex colors.Note that one can reduce the König representation by omitting the vertices which represent the sites which are only contained in one bond within the cluster [50].The number of the omitted sites are easily deduced from the bond-types.We do not distinguish flux and charge-sites within the hypergraph representation, because they are treated on equal footing in the operators τ n,b (α) and also in the unperturbed Hamiltonian.So basically, the problem is decomposed into stabilizer eigenvalues which are relevant for the energetics, and spin-background variables, which we introduce to obtain the correct phases.It is clear that the additional phase depends only on the action of the Pauli matrices on the reference state.This is correctly taken into account, once we know which Pauli matrices act on which spin-background variable and in which order.From the edge-colors in the König representation it is clear, which sites correspond to the spin-background.The order is given by the sequence of operators τ n,b (α) .Regarding the quasiparticles, the graph representation might result in graphs for which it is not clear whether vertices correspond to charge or plaquette sites.However, energetically this does not matter.In summary, we only separate the energetic degrees of freedom from the spin-background within the graph representation exploiting that elementary fluxes and charges have the same energy.
Another important question is how obtain the embedding factor N (E, C).Typically, the ground-state energy is calculated per site, so N (E, C) is the number of sub-clusters in E per spin-site.Here we calculate it per spin-site, i.e., modulo translational symmetry.To this end we count all connected sub-clusters of the system which contain a given spin-site s i and belong to E. Clearly, within this set there are still sub-clusters related by non-trivial translations.However, it is also clear that the number of translations of a sub-cluster which contain the spin-site s i id equal to the number of spin-sites within the sub-cluster.So in order to obtain N (E, C), one can just count the number of sub-clusters in E which contain s i , and divide by the number of sites of any sub-cluster in E. One can also require that for any two sub-clusters in E there exists an isomorphism taking the vertices representing the site s i to each other.Note that we still require s i to be on the sub-clusters.To obtain N (E, C) the number of these sub-clusters is divided by the size of the orbit of the spin-site s i under the action of the automorphism group of the representing graph.

Single-particle excitation gap
For the single particle irreducible matrix elements only minor modifications are necessary in the graph representation.Recall that we are aiming to calculate irreducible matrix elements of the form with ⃗ δ = ⃗ r 2 − ⃗ r 1 .In Subsect.4.2 we established that using appropriate conventions for the calculations we can execute the calculations on finite clusters, and they are not affected by any degrees of freedom outside the cluster.So we can encode the product state on the cluster using vertex colors for all vertices which represent sites.Note that these vertices must still be distinguishable from the vertices which represent bonds.Depending on the local state at a where the states label the values of x , z or si at the respective site.A white filling of the vertex means, that the eigenvalue associated with the site is zero within both states |⃗ r 1 〉 and |⃗ r 2 〉, while an orange vertex indicates that the respective eigenvalue changes from 1 to 0. A specific example for this representation is shown in Fig. 8.In summary we distinguish seven different types of vertices with different colors: Three different types of vertices representing bonds b (x) , b ( y) , b (z) and four different types of site vertices representing changes of the local states as indicated in Eq. (58).Next, we adress the question, which sub-clusters we need to consider.The obvious answer is every sub-cluster c with 〈⃗ First of all, it is clear that a sub-cluster can only contribute if it contains all stabilizer sites, which change their local state within the matrix element.Next, the sub-cluster should host excitations or the initial state |⃗ r 1 〉 should have a non-trivial spin-background.Interestingly the spin-background of |⃗ r 2 〉 does not matter.In order to determine the state |f〉 resulting from the action of a sequence of operators τ n,b (α) on a state |⃗ r 1 〉 the spin-background of the state |⃗ r 1 〉 is crucial.Using an appropriate reference state the product of operators τ n,b (α) yields the result in terms of a phase factor [i n a ], and a state equivalent to a canonical state |⃗ r ′′ 〉.So the spin-background of the state |⃗ r 2 〉 is not important, as This also implies that the spin-background of |⃗ r 2 〉 does not play a role in the representation and the selection of the sub-clusters.As a consequence, we can ignore the spin-background of |⃗ r 2 〉 for the graph representation, while we still need to include the spin-background of |⃗ r 1 〉 With this the König representation from Fig. 8 is replaced by the representation in Fig. 9.
For the stabilizer degrees of freedom we keep the convention defined in Eq. ( 58) so we still distinguish seven different types of vertices.This way of coloring the vertices has the advantage that the rules which clusters need to be taken into account are formulated more easily.For offdiagonal elements 〈⃗ r 2 , | H eff − E 0 |⃗ r 1 , 〉 with ⃗ r 2 ̸ = ⃗ r 1 all sites represented by green or orange vertices need to be included in a cluster, otherwise its reduced contribution vanishes.Instead, for the diagonal element the reduced contribution of a cluster is zero, if no site represented by a black vertex is contained in the cluster.A specific cluster on which the winding process of a flux around a charge in Fig. 5 occurs is represented by the graph in Fig. 10.While this graph does not contain any vertices corresponding to excitations the non-trivial spin-background suffices to take the effects of the charge (which is not on the cluster) into account.Note that this graph can only play a role within matrix elements 〈⃗ r Note that also clusters which do not enclose a charge can have a non-trivial spin-background in such a hopping element.Currently, we simply consider the respective equivalence classes within the calculations.It is interesting whether a scheme to sort out these clusters can be devised and how it affects the efficiency of the method.

Series for ground-state energy and excitation gaps
Using the described hypergraph decomposition, we calculated the ground-state energy per spin up to order 10 using the pCUT method.The ground-state energy, which is symmetric under the exchange of h x and h z due to self duality, can be written down by summarizing certain terms using the notation S j = h The one-particle excitation energies have been calculated in order 9 for arbitrary field directions and in order 10 for a parallel field direction with h y = 0.Here we only focus on the charge and flux gap which is always located at ⃗ k = (0, 0).Due to self-duality, we can constrain ourselves to the charge gap ∆ which reads for arbitrary field direction Due to reduced computational time, in the parallel field case with h y = 0, the charge gap was calculated up to order 10 The flux gap ∆ can be again obtained by exchanging h x and h z in these expressions.These results agree with previous high order series expansions [25,53].In both the arbitrary and parallel field cases we were able to exceed the previously known results by one order.

Conclusions
In this work we presented how to treat non-local anyonic statistics with a high-order linkedcluster expansions using a full hypergraph decomposition.This we exemplified for the topological phase of the toric code in the presence of a uniform field possessing Abelian charge and flux quasi-particles.Technically, we used the pCUT method to obtain series expansions for the ground-state energy and the one-quasiparticle charge and flux gap.Note that also other perturbative techniques like Takahashi's perturbation theory [72,73] or Löwdin's partitioning technique [74,75] can be applied to calculate these properties of the perturbed topological phase.
The improvement in terms of the achieved maximal perturbative order compare to calculations done with Entings finite lattice method [21,22,25] is only moderate.Nevertheless, the explicit treatment of the non-locality of the fractional statistics of charges and fluxes within a full hypergraph decomposition is an asset on its own.Further, there will be a benefit of a full hypergraph decomposition when performing similar calculations for three-dimensional perturbed topological phases like the ones in the three-dimensional toric code [38][39][40] or the fractonic X-Cube model [76] in the presence of external perturbations.Finally, while we expect that extensions to other Abelian topological phases can also be treated in the described framework, extensions to perturbed non-Abelian topological phases pose an interesting and open challenge.

A Identifying non-contributing clusters or isomorphism classes
The selection rules given in [57,63] can be adapted for the problem presented in this paper.We will shortly describe how this can be done but refer to [50] for a longer discussion in the context of hypergraph expansions.We will explain the heuristics for the vacuum.A generalization for hopping elements is straightforward.

A.1 Vertex Degree
We are considering the degrees of the vertices which represent stabilizer operators.Any such vertex which has an odd degree within the full König representation induces that the perturbation has to act more than once on an adjacent bond.Of course such a vertex can share this bond with other vertices of odd degree, though the maximum number of sites in a bond c max poses an upper limit to this.Here c max = 4 for the general field and c max = 2 for a parallel field (h y = 0).So having m vertices of odd degree which represent a stabilizer eigenvalue we need to act at least ceil(m/c max ) additional times with the perturbation, where ceil denotes rounding up to the next integer.For a cluster C with n b bonds and m such vertices any process contributing to κ(C) has to act at least n b + ceil(m/c max ).So at perturbation order k the contribution κ(C) of a cluster C will be zero if It is obvious, that if a cluster fulfills (65) all clusters created adding further bonds will fulfill the inequality as well.This procedure is applicable to clusters and isomorphism classes as well.It is very powerful because it prunes whole subtrees during the search for new clusters.On the cluster level this scheme can be extended taking into account the actual bonds of the system.This can be achieved by storing whether two sites are contained in a at least one of the bonds simultaneously for any pair of sites.Note that here we are considering all the bonds in the entire system.Once we have the vertices of odd degree which represent the stabilizer eigenvalues, we know which of these vertices have to be in different bonds.As an example suppose we have four such vertices in a cluster C, and all of these vertices have to be in different bonds.With this we mean, that there is no bond containing two of these vertices in the entire system (including also the bonds in C).Then from the simple estimate above we conclude that the cluster contributes in order n b + 1 (considering h y ̸ = 0), where again n b is the number of bonds in the cluster C. Instead, with the refined estimate we conclude that it contributes only in order n b + 4. How one can estimate a lower bound for the number of different bonds needed using graph coloring algorithms is explained in [50].Let us conclude with the remark that also this estimate can be used to prune entire branches of the search tree.It is advisable to use both of these techniques during the generation of the clusters, in order to improve the performance of the cluster generation.

A.2 Doubling Edges
Another possibility to sort out non-contributing isomorphism classes is a technique presented by He et al. [63].We consider a cluster C representing a given isomorphism class and try to obtain a cluster where all vertex degrees (again here we consider only vertices which represent stabilizer degrees of freedom) are even by duplicating bonds within the cluster C. The smallest number of bonds of such a cluster is the order in which the isomorphism class can contribute.We suggest applying this technique after the cluster generation on the level of the isomorphism classes.This reduces the number of considered equivalence classes further and avoids superfluous calculations.

1 Figure 1 :
Figure1: Left: The spin-1/2 degrees of freedom are depicted in light blue at the edges of the square lattice.The operators X and Z are illustrated in red or blue respectively.Blue (red) circles indicate that the illustrated operator acts with σ z (σ x ) on the respective spin.Right: Two examples for string operators of the type S x (S z ) in red (blue).Actions of σ x (σ z ) on spin degrees of freedom are represented by red (blue) circles.The flux (charge) excitations at the end of the string operator S x (S z ) are illustrated as blue (red) squares labeled with −1 indicating the negative eigenvalue z (x ) of Z (X ).

Figure 2 :
Figure2: Canonical one-particle states.Left: A single charge excitation is obtained by creating two charges from the vacuum and moving one to infinity.This procedure results in a semi-infinite string attached to the charge (indicated in blue).Right: A single flux excitation is obtained by creating two fluxes from the vacuum and moving one to infinity.This procedure results in a semi-infinite string attached to the flux (indicated in red).

Figure 3 :
Figure 3: Left figure illustrates a non-canonical one-charge state.The product of a contractible loop operator of σ z with a non-canonical string operator (left) is equal to a canonical string operator (right).Recall that a contractible loop of σ z is equivalent to a product of operators Z .

Figure 4 :
Figure 4: An illustration of the three different bond types b (α) i with α ∈ {x, y, z}.These bonds contain the respective charge and flux sites as well as the spin s i which are affected by the action of Pauli matrices σ x i , σ z i and σ y i .

Figure 6 :
Figure 6: The König representation of the different bonds with the contained sites is shown.The double circles represent the bonds, and the single circles represent the sites, namely the star and the plaquette eigenvalues and the spin background variable.Note the site corresponding to the spin-background variable is connected to the double circle by a thick red edge.

Figure 7 :
Figure 7: Left: Example for the König representation of a cluster.The bonds are represented by double circles, where the bond types are explicitly given.The different labels x, y, z as well as double and single circles correspond to different vertex colors.The vertices corresponding to spin-sites are explicitly distinguished as they are connected by red edges to the vertices representing the bonds.Right: In the reduced König reprsentation the vertices from the vertex-part are omitted if they have degree one, reducing the complexity of the graphs representation.

9 Figure 8 :
Figure 8: Left: Illustration of the König representation for the hopping element 〈10000000; 11| H ceff − E 0 (c) |01000000; 01〉.The i-th number in the states refers to the vertex with number i in the König representation.Note that the enumeration of site vertices is not part of the actual representation, whereas the label of the bondvertices is.We do not explicitly distinguish charge and spin sites, e.g., we cannot deduce which types of sites the vertices 5, 6, 7 represent.The label x, y, z of the bond vertices is part of the representation.Right: By omitting the white site-vertices of degree one we again obtain a reduced representation.

9 Figure 9 :
Figure 9: Illustration of the König representation of a hopping element, without taking into account the spin-background of the state |⃗ r 2 〉.In contrast to Fig. 8 the vertex 8 is now black instead of orange.

Figure 10 :
Figure 10: A graph which can only contribute to diagonal matrix elements 〈⃗ r, | H eff − E 0 |⃗ r, 〉.The spin-background distinguishes this graph, representing a cluster which encloses a charge from a graph which represents a cluster which does not enclose a charge.