Topological spin liquids: Robustness under perturbations

Mohsin Iqbal,1,2 Helena Casademunt ,3,4 and Norbert Schuch1,2 1Max Planck Institute of Quantum Optics, Hans-Kopfermann-Str. 1, D-85748 Garching, Germany 2Munich Center for Quantum Science and Technology, Schellingstr. 4, D-80799 München, Germany 3Department of Physics, Princeton University, Princeton, New Jersey 08544, USA 4Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA


I. INTRODUCTION
Topological spin liquids (TSL) are exotic phases of matter where a system does not order magnetically despite strong antiferromagnetic interactions, but rather topologically, i.e., in its global entanglement. The interest in these systems stems from their unconventional properties, such as anyonic excitations with fractional charge and nontrivial statistics, and a ground space protected by its global entanglement [1][2][3]. However, TSLs are notoriously difficult to identify, both in theory and in experiment, as candidate systems often exhibit close competition between a number of different phases. In order to robustly realize these phases, it is therefore essential to understand how sensitively they react to perturbations which can induce a breakdown of topological order.
The paradigmatic example of a spin liquid is the resonating valence bond (RVB) wave function on the kagome lattice, which consists of a "resonating" superposition of all possible ways to cover the lattice with nearest-neighbor singlets [4,5]. It forms a physically motivated low-energy ansatz for Heisenberg-type models, and appears as the exact ground state of a local model with topological order [6,7]. However, the nonorthogonality of different singlet configurations makes RVB models hard to analyze. To mitigate this difficulty, dimer models have been studied instead, where different singlet configurations are taken to be orthogonal [8]. The resulting kagome dimer model is an RG fixed point and thus significantly easier to analyze [9]. However, it is unclear to what extent results derived for the dimer model still apply to the RVB state, where the nonorthogonality of singlets induces a finite correlation length, making it doubtful whether the robustness of the RVB state can be understood from studies performed on the dimer model.
In this paper, we study and compare how sensitively the RVB spin liquid and the quantum dimer model react to noise, and which level of perturbations they can tolerate before their topological order breaks down. We consider both magnetic fields and lattice anisotropies, corresponding to doping with the two elementary topological excitations: spinons and visons. We find, rather surprisingly, that in both cases the RVB model exhibits essentially the same stability as the dimer model despite its nonzero correlation length. This suggests that the dimer model is more accurate in modeling spin liquid physics under perturbations than one might have naively assumed.
To understand the mechanism behind this unexpected result and its range of applicability, we microscopically analyze the structure of anyon correlations using tensor networks. Diverging anyon correlations indicate a closing gap, driving a phase transition through anyon condensation; the finite spinon correlations in the RVB would thus indeed suggest a decreased robustness. However, as our analysis reveals, there is a universal mechanism underlying the surprising robustness of the RVB model: It arises from symmetries which protect specific correlations, and is thus independent of the specific perturbation but rather a universal feature of the RVB spin liquid.  1. (a) Dimer pattern (blue) and its arrow representation (red); arrows point into the triangle with the dimer. (b) The difference between any allowed arrow pattern (green) and the reference pattern (red) is in one-to-one correspondence to loop patterns. (c) Construction of "dual tension" doping; cf. text.
Concretely, we find that the nonorthogonality of singlets induces dominant spinon correlations without a separate vison correlation scale; since vison doping only increases the latter, the response is unaffected by the spinon length scale. The reason for the robustness to spinon doping is more subtle; we assess it through a combination of analytical and numerical study. It reveals that the spinon correlations exhibit a twofold degeneracy in addition to the spin doublet. It originates in the two different ways to construct spinon correlations-either through the overlap of two states with one spinon each, separated by some distance , or through the overlap of the original doped RVB with a state with two additional spinons placed at a separation . We show that the reflection symmetry of the RVB rules out any correlation between spinons of opposite spin, also at finite doping. This, together with the symmetry of these overlaps under ket ↔ bra exchange, implies that the fourfold multiplet splits under doping into sectors labeled by (i) their spin ± 1 2 and (ii) a fixed relative phase ±1 between the spinon at either position being in the ket or bra vector in the overlap, respectively. As we show, this ±1 phase label is protected by the reflection symmetry of the lattice and thus stable to perturbations which respect that symmetry. We find that the correlations which are initially enhanced by magnetic fields and those which drive the phase transition live in different sectors which are protected by the lattice symmetry; that explains why the presence of initial spinon correlations in the RVB state has no effect on its robustness to magnetic fields.

II. RVB AND DIMER MODEL
The RVB state is constructed as follows. First, we define a dimer covering as a full covering of the lattice with pairs of adjacent vertices, termed dimers [blue in Fig. 1(a)]; we denote dimer coverings by D, and the set of all coverings (with PBC) by D. Next, replace each dimer by a singlet |σ = 1 √ 2 (|↑↓ − |↓↑ ) (with counterclockwise orientation around triangles), we call the resulting state |σ(D) . The RVB wave function is then |RVB = D∈D |σ(D) . In studying the physics of RVB wave functions, so-called dimer models |dimer = D∈D |D are frequently used, where the {|D } D∈D define an orthonormal basis [8]. Replacing singlet configurations by orthogonal dimer configurations makes these models easier to analyze, but can also affect their physics. One way to explicitly construct a dimer representation is to start from the RVB wave function and attach arrows to each vertex [10] which for any dimer configuration D point into the triangle where the adjacent dimer lies [ Fig. 1(a)]. These arrows can be treated as quantum degrees of freedom or "arrow qubits" with basis states {|a ↑ , |a ↓ } (arrow pointing into either of the two adjacent triangles); denoting the corresponding global arrow configuration by |A(D) , we obtain a local representation |dimer = D∈D |σ(D) |A(D) of the dimer model. One advantage of this representation is that it allows one to continuously interpolate between the dimer model and the RVB state, by choosing a nonorthogonal arrow basis |a ↑/↓ = (1 ± λ)|0 + (1 ∓ λ)|1 and tuning λ ∈ [0; 1]. It can be proven that along the whole interpolation, the system has a parent Hamiltonian with a fourfold degenerate ground space with topological features, and numerical study shows that the correlation length along the interpolation stays finite, placing both models in the same (topological) phase without any conventional order [6].
The dimer model can be proven to be a topological fixed point model using the arrow representation introduced above. First, given any classical configuration |A(D) , we can disentangle the singlets |σ(D) by local unitaries (conditioned on the adjacent arrow qubits) and bring them to a fiducial state, leaving us with a superposition A∈A |A of all allowed arrow configurations A ≡ A(D); these are precisely those with an even number of inpointing arrows (a Z 2 constraint) [10]. By fixing a "reference configuration" A 0 of arrows [and thus a reference configuration D 0 of dimers; Fig. 1(a)], every arrow configuration A ∈ A is characterized by those vertices where the arrows in A differ from A 0 . These vertices satisfy a Z 2 Gauss law on each triangle and thus describe closed loops L on the dual honeycomb lattice [ Fig. 1(b)]; this establishes a one-to-one correspondence between loop configurations L and arrow configurations A. The dimer model is thus locally equivalent to an equal-weight superposition of all loop configurations on the dual honeycomb lattice, which is nothing but the topological Z 2 toric code model [11]. In fact, we can think of the dimer and RVB model as being constructed from a loop model to which we apply a sequence of local operations which replace the loops by arrows, add singlets as prescribed by the arrows, and finally (partially) erase the arrow pattern by applying to each arrow qubit; note that while the first two steps are local unitaries/isometries, the last step E θ is a nonunitary ("filtering") operation which induces a finite correlation length and in the limit λ → 0 becomes singular.

III. DOPING THE RVB WAVE FUNCTION
Subjecting the RVB or dimer model to external fields induces doping with elementary excitations: spinons or visons. Spinons are obtained by breaking up a singlet (or dimer) and replacing it by two separate spins (w.l.o.g., two up-spins), which can subsequently separate due to a kinetic term. Visons, on the other hand, correspond to a local disbalance in the relative weight of different singlet (or dimer) configurations (equivalently, loop configurations).

115101-2
In order to study how a finite density of excitations affects the topological order in the RVB or dimer model, we extend the ansatz to include a tunable quasiparticle doping. Let us start with vison doping: Here, we select a reference dimer pattern D 0 [=arrow configuration A 0 , Fig. 1(a)], and adiabatically increase the relative weight of the reference configuration through a filtering F θ v = 1 − θ v |r r| (with |r the opposite of the reference state) applied to each arrow qubit before the application of E λ . This directly translates to a "string tension" in the underlying loop model (defined relative to D 0 ) which suppresses longer loops and gives rise to doping with magnetic (vison) excitations; for the dimer point λ = 1, the two doping models are unitarily equivalent.
For spinon doping, we introduce two ansatzes. The firsttermed "dual tension"-maps to electric excitations (broken loops). This corresponds to flipped arrows, obtained by applying G θ s = 1 + θ s σ x Two inpointing arrows are mapped to a singlet, and one (three) to |↑ (|↑ ⊗3 ) [ Fig. 1(c)]. However, the resulting ansatz does not correctly reproduce the effect of a local field in lowest order-breaking a singlet into a pair of spinons-as it also yields four-spinon terms. We therefore introduce a second ansatz ("spinon pairs") by replacing the singlets in |σ(D) with a pair of spinons |↑↑ , i.e. changing each singlet to 1 Unlike the other ansatz, this correctly captures the expected behavior in leading order. We have also found that it performs significantly better as a variational ansatz for the Heisenberg model with magnetic field. We therefore focus on it, and discuss the other ansatz in Appendix B.

IV. RESULTS AND ANALYSIS
We will now study the response of the RVB spin liquid to different fields. For all simulations, we have expressed the wave functions as projected entangled pair states (PEPS) [6,[12][13][14]; see Appendix A. We use standard numerical PEPS methods (boundary MPS [15][16][17]) which allow us to evaluate physical observables in the thermodynamic limit, as well as to extract a correlation length from the boundary MPS. Using the techniques in Refs. [17,18], we can moreover extract correlation lengths for each anyon sector, which label the decay of correlations between a pair of anyons of a certain type. This allows us to microscopically analyze how the system is driven into a trivial phase due to doping with some anyon a, causing it to condense-in this process, the mass gap of a decreases until it eventually vanishes and it becomes favorable to have a macroscopic number of a anyons in the ground state: Anyon a has become condensed, leading to a breakdown of topological order [19]. In order to probe the anyon mass m a , we can study the anyon-anyon correlation length ξ a ∼ 1/m a ; a diverging ξ a thus indicates condensation of anyon a.
We start by considering doping with visons due to lattice anisotropies which drive the system into a vison condensed phase [i.e., a valence bond crystal (VBC)] [20]. θ v . We find that the critical point θ crit v (λ) is essentially independent of the interpolation λ between the dimer and the RVB model. Figure 2(b) shows the response to a doping θ v , defined as the difference between the Heisenberg energies on inequivalent edges, confirming that the RVB and the dimer model behave essentially the same way. We can understand this behavior by considering the correlations (mass gaps) for the different anyon types in the RVB state as a function of θ v [Fig. 2(c)]. We find that at the RVB point, the dominant length scale is given by spinon-spinon correlations ξ s [21]. Vison correlations ξ v , while present, are only on the order of the topologically trivial correlations ξ t ≈ ξ v , which in turn are roughly ξ t ≈ ξ s /2, and thus understood as arising from correlations between two pairs of spinons (which are topologically trivial). That is, the dominant length scale in the system arises from spinons, while visons do not exhibit independent correlations on their own. As we increase the doping θ v , genuine vison correlations start building up, but the overall correlation length remains dominated by the spinon correlations, which do not respond to the vison doping; only very close to the phase transition (θ v ≈ 0.18), the vison correlations start to exceed the spinon correlations and diverge at the phase transition. Next, we consider the effect of magnetic fields, amounting to doping with spinons (the "spinon pairs" ansatz). Figure 2(d) shows the phase diagram as a function of the doping θ s and the dimer-RVB interpolation λ. Surprisingly, we find that despite the dominant spinon correlations in the RVB state, the phase boundary stays almost constant (the RVB is even slightly more robust). Again, studying the response m z = 1 N σ z vs θ s , Fig. 2(e), we find two very similar curves, and the RVB shows a smaller response in the relevant regime; as expected, the phase transitions coincide with maximal susceptibility. 1 These findings are rather counterintuitive, given the dominant spinon correlations in the RVB. To analyze this, we consider the correlations by anyon type [Fig. 2(f)]: We find that the spinon correlation dominates throughout, but it decreases after an initial increase, and exhibits a sharp kink around θ s ≈ 0.99 after which it diverges. To study this further, we consider the full spectrum of correlations in Fig. 2(g), where we make two noteworthy observations. First, the leading spinon correlation is fourfold degenerate, where one would have naively expected a spin- 1 2 doublet. Moreover, the correlation in this quadruplet which dominates at small doping is different from the one which finally drives the phase transition-the two lines exhibit a sharp crossing at θ s ≈ 0.99, suggesting they are distinguished by some symmetry. Indeed, such a symmetry protection could explain the surprising robustness of the RVB model, since the correlations in the sector driving the phase transition initially decrease under doping.
To analyze this further, we first consider the origin of spinon correlations. Correlation functions are overlaps ψ |ψ , where both |ψ and ψ | are RVB states possibly doped with spinons-that is, they are a sum of all singlet coverings, except for one or two static locations where |↑ spinons are placed. We will call |ψ the "ket layer" and ψ | the "bra layer." Spinon-spinon correlations are obtained by summing over all overlaps of singlet patterns with the two spinons fixed (Fig. 3), where each pattern yields an amplitude determined by the loop lengths, and the sign follows from the singlet orientations. There are two types of such correlations: 1 Note that m z does not directly measure the spinon density for the RVB, since nonzero contributions also arise from pairs of spinons connected by singlets in ψ|σ z |ψ .
A ↑ ket spinon (i.e., a single up-spin in the ket layer) can be correlated either to a ↑ bra spinon [ Fig. 3(a)], or a ↓ ket spinon [ Fig. 3(b)] (since S ket z,total = S bra z,total ), and correspondingly for a ↓ ket spinon. We can now explain the fourfold degeneracy: In the S total z = + 1 2 sector, either a ↑ ket spinon is correlated with α↑ bra + β↓ ket in a fixed superposition, or independently ↓ bra with a superposition α↓ ket + β↑ bra . This yields a twofold degeneracy, while another factor of 2 arises from the spin-1 2 doublet.
However, there is more structure to the spinon-spinon correlations: Every overlap pattern in Figs. 3(a) and 3(b) has a "twin" under reflection about the indicated vertical axis. For the singlet orientation chosen, all singlets are flipped under reflection: Thus, overlaps for odd length paths-which connect ket-ket and bra-bra spinons-change their sign under reflection and thus cancel: 2 In the RVB state, only spinons with same S z (equivalently, in opposite layers) can be correlated. Remarkably, this property is preserved under doping [ Fig. 3(b)]: Since spinon pairs are symmetric under reflection and do not flip the spin, the paths containing an odd (even) number of spinon pairs which don't cancel with their reflections are exactly those which correlate spinons with the same S z in the same (opposite) layer. Together with the ket↔bra symmetry, this implies that the spinon-spinon correlations can be labeled by sectors ↑ ± := ↑ ket ± ↑ bra and ↓ ± := ↓ ket ± ↓ bra (that is, sectors which have overlap only with the corresponding superposition of spinons). At the RVB point, they all appear with equal amplitude but opposite phase such that the equal-layer correlations cancel. Crucially, as shown above, this symmetry label is protected by reflection symmetry even at finite doping, preventing mixing. Our analysis is confirmed by numerically computing the matrix elements of ± with the different correlation sectors. 3 Moreover, analysis of the data shows that the correlations which increase initially are in the − sectors, while those in the + sectors decrease, whereas the phase transition is driven by a diverging correlation length in the ↑ + sector. However, coupling between the two sectors is prohibited by reflection symmetry. This explains why the system responds to spinon doping with an increased correlation length, and yet, this does not come with a decrease in robustness of the topological phase: The two phenomena take place in different symmetry sectors. A qualitatively similar behavior is observed for the "dual tension" doping: The spinon correlation spectrum again splits in the ± basis, and we find only a very weak effect on the robustness.
While our analysis is based on a specific doping model, this is in fact a universal behavior, since any perturbation which results in a breaking of singlets into pairs of spinons will have the same effect to leading order. This implies that the initial splitting will again be in the ± basis, with the − correlations being dominant for small doping, while the ↑ + correlations drive the phase transition. As long as the perturbation respects the lattice symmetry, these correlation sectors cannot mix, and we therefore expect a qualitatively similar behavior for a general perturbation which induces doping with spinons.

V. CONCLUSIONS
We have studied the robustness of the RVB and the dimer model on the kagome lattice to perturbations using PEPS. We have found that despite the nonorthogonality of different singlet configurations, the RVB spin liquid exhibits the same robustness to perturbations as the dimer model. For lattice anisotropies (doping with visons), we traced this back to the fact that the length scale induced by the nonorthogonality of singlets gives rise to spinon correlations but does not directly affect the physics of visons. For magnetic fields (doping with spinons), we showed that the robustness arises from a protection of the RVB state due to the reflection symmetry of the lattice, which separates the initially dominating spinon branch from the branch which ultimately drives the phase transition. Our results reveal a surprising universal robustness of the RVB spin liquid against perturbations, highlighting its role as a candidate for the realization of a stable gapped spin liquid. Here we describe how the different doping mechanisms introduced in the paper can be described as tensor network states, also termed projected entangled pair states (PEPS).

String tension and dual tension
Let us first describe the PEPS for the RVB with vison and spinon dopings constructed from doping of the underlying loop (or arrow) model, i.e., string tension or dual string tension. This construction will consist of three layers, stacked on top of each other. The lower layer is a PEPS for the loop model with the corresponding (dual) tension. On top of that layer, we apply a projected entangled pair operator (PEPO) which transforms this loop model into the corresponding dimer model. Finally, in a last step we apply local filtering operations (as introduced in the main text) to the arrow degrees of freedom in the dimer model, which allows one to interpolate to the corresponding RVB state.
The PEPS for the first layer-the loop model with tension-consists of two types of tensors: one vertex tensor (without physical legs) and an on-site tensor (which carries the physical index). The vertex tensor has three legs, each of dimension two. We use a computational basis to express the presence or absence of a loop string on the link, and due to the Z 2 constraint, the vertex tensor is restricted to four nonzero entries, on physical legs. (In the main text, we restrict to the cases where one of the θ • ≡ 0). The second layer is a PEPO which maps loop configurations on the honeycomb lattice (including broken loops in the case of dual tension) to dimer configurations on the kagome lattice (including monomers in the case of broken loops). It again consists of two types of tensors. The first is a triangular tensor without physical index, which has three indices of dimension three each: where the tensor σ = diag( 1 √ 2 ( 0 −1 1 0 ), 0) models an oriented singlet and the last four terms correspond to different spinon configurations. The second tensor acts on each site: It takes the loop configuration as an input in index a, and outputs a physical spin p and an arrow index d, and is built such as to pick the physical qubit from either of the two adjacent virtual indices [and thus triangular tensors (A4)], as prescribed by the reference configuration: By design, this tensor is not symmetric in the virtual indices i and j, and we use an arrow pointing to the index j to label its orientation. The PEPO [ Fig. 4(b)] is now obtained by assembling (A4) and (A5) in a hexagonal structure (yielding spins on a kagome lattice), where the arrows need to be oriented such that setting a ≡ 0 yields the reference configuration.
The tensor network for the doped dimer model is now obtained by stacking the two tensor networks in Fig. 4(a) (for the honeycomb loop model) and Fig. 4(b) (for replacing loops with dimers), where the gray indices (labeled a) are contracted. The resulting tensor network allows one to tune the doping with spinons and visons by changing the parameters θ s and θ v in the deformation tensor (A3).
Finally, by applying a filtering, on the arrow qubits d, we can continuously interpolate between the dimer and RVB models also with doping.

Doping with spinon pairs
The tensor network to implement the "spinon pair" doping is obtained by modifying the second layer of the preceding construction. There is no longer a need for the first layer (the loop model), since the loop constraint is already contained in the dimer model (see also the original construction for the RVB and dimer PEPS [6]). First, the tensor (A4) is modified such that each index has dimension four: Here, the tensor σ = diag( 1 √ 2 ( 0 −1 1 0 ), 0, 1) encodes either a singlet (in the first two basis states) or the presence of a spinon pair (in the new fourth degree of freedom). Correspondingly, the on-site tensor is also changed to project to the spinon degree of freedom with a tunable weight of θ s , accompanied by a third state d = 2 of the arrow qubit (the basis state |a s ): where the arrows are oriented as before. In order to interpolate to the RVB state, we can erase the information on the dimer Comparison of Heisenberg energy vs magnetization for the two models. This illustrates that the two models already differ in the perturbative regime (close to the origin); only the "spinon pairs" ansatz correctly captures the physics of breaking a singlet into a nearest-neighbor pair of spinons in leading order. (c) Correlation functions by spinon sector for the "dual tension" spinon doping; we observe the same characteristic features as in Fig. 2(f), which are protected by the same symmetries (in particular lattice reflection) as discussed in the main text. (d) The phase diagram of the dimer-RVB interpolation θ s for the "dual tension" doping λ again exhibits only a weak dependence of the phase transition on changing orthogonal dimers to nonorthogonal singlets. indices by applying a deformation, The final tensor network is identical to Fig. 4(b), but without the gray indices.

APPENDIX B: SPINON DOPING WITH "DUAL TENSION"
In this Appendix, we report the results for the model with spinon doping constructed through dual string tension. Figure 5(a) provides a comparison of the variational energy for the "dual tension" and the "spinon pairs" ansatz for the Heisenberg Hamiltonian with a transverse field, (with eigenvalues S z i = ± 1 2 ). We find that the energy for the "spinon pairs" ansatz is significantly lower, providing a first reason why we chose to consider it as our primary ansatz for spinon doping. Figure 5(b) provides further insight into this. It shows the relation between Heisenberg energy S i · S j and magnetization S z i : We see that for the same magnetization, 115101-6 the "dual tension" ansatz has a significantly higher Heisenberg energy, which qualitatively means that it requires to break up a correspondingly larger number of singlets to achieve the same magnetization. This effect can be clearly observed in the perturbative regime, i.e., small magnetizations, where the "dual tension" ansatz has a significantly higher slope. We can understand this effect qualitatively in a semiclassical picture: Breaking up a singlet into a pair of spinons leads to a change E = 1 in Heisenberg energy, since one singlet is replaced by a |↑↑ state. On the other hand, the scenario in the "dual tension" construction where flipping an arrow yields four spinons (three on a triangle, and one adjacent vertex) gives rise to a total of four Heisenberg terms having an energy + 1 4 , while before, half of them had energy 0 and half − 3 4 , implying E = 2.5 (or 1.25 per pair of spinons).
Despite these differences, the study of correlations by anyon sectors, Fig. 5(c), yields a qualitatively very similar behavior to the case of "spinon pair" doping. In particular, we again observe an additional twofold degeneracy in the spinon sector (on top of the spin- 1 2 multiplet) which is protected by the lattice symmetry, and which separates the correlations responsible for the phase transitions from those initially responding to the doping; this highlights the fact that this, as we have shown, is a universal effect. Yet again, this symmetry protection is reflected in a rather weak dependence of the phase transition on interpolating between the doped orthogonal dimer model and the doped RVB state [ Fig. 5(d)].