Gapped $Z_2$ spin liquid in the breathing kagome Heisenberg antiferromagnet

We investigate the spin-1/2 Heisenberg antiferromagnet on the kagome lattice with breathing anisotropy (i.e. with weak and strong triangular units), constructing an improved simplex Resonating Valence Bond (RVB) ansatz by successive applications (up to three times) of local quantum gates which implement a filtering operation on the bare nearest-neighbor RVB state. The resulting Projected Entangled Pair State involves a small number of variational parameters (only one at each level of application) and preserves full lattice and spin-rotation symmetries. Despite its simple analytic form, the simplex RVB provides very good variational energies at strong and even intermediate breathing anisotropy. We show that it carries $Z_2$ topological order which does not fade away under the first few applications of the quantum gates, suggesting that the RVB topological spin liquid becomes a competing ground state candidate for the kagome antiferromagnet at large breathing anisotropy.


I. INTRODUCTION
The Resonating Valence Bond (RVB) state, defined as an equal weight superposition of (non-orthogonal) nearest neighbor (NN) singlet (or dimer) coverings, was first proposed by Anderson [1] to describe a possible spin liquid ground state (GS) of the S = 1/2 antiferromagnetic Heisenberg (HAF) model on the triangular lattice. Later on, it was also introduced as the parent Mott state of high-T c superconductors [2]. Several works [3][4][5][6] have demonstrated that NN RVB states defined on triangular and kagome lattices are gapped spin liquid states with Z 2 topological order, and GSs of local parent Hamiltonians [4,7].
Spin liquid behaviors are expected in two-dimensional (2D) frustrated quantum magnets where magnetic frustration prohibits magnetic ordering at zero temperature. The spin-1/2 Heisenberg antiferromagnet on the kagome lattice (KHAFM) is believed to be the simplest archetypical model hosting a spin liquid GS with no Landau-Ginzburg spontaneous symmetry breaking. However, the precise nature of this spin liquid is still actively debated. While the HLSM theorem [8] excludes a unique GS separated from the first excitations by a finite gap (so-called "trivial" spin liquid), a gapless spin liquid [9][10][11][12] or a gapful topological spin liquid (of the RVB type) [13][14][15] are the two favored candidates.
An important aspect is to understand the stability of the spin liquid GS against various perturbations, such as long-range interactions or different anisotropies. Beyond being of interest by itself, it might also yield alternative ways to assess the nature of the ground state of the isotropic KHAFM by allowing to adiabatically connect it to a limiting case which might be easier to study. An important case of such perturbations is the HAF model on the kagome lattice with anisotropy, which can be written as with 0 ≤ γ ≤ 1 (where S z = ± 1 2 ). The Hamiltonian H(γ), except at γ = 1, explicitly breaks the inversion symmetry between the strong (or right-pointing) and the weak (or left-pointing) triangles of the kagome lattice (Fig. 1). The anisotropic model (1) (also referred to as "breathing" HAF [16]) has gained additional relevance because recent studies have shown a realization of (1) for particular values of γ in a vanadium-based compound [17][18][19]. Moreover, in the limit of strong anisotropy, γ → 0, it can be mapped to a simpler model with two spin- 1 2 degrees of freedom per site, similar to a Kugel-Khomskii model [20].
The Hamiltonian (1) has been studied using different numerical methods: In Ref. 21, Gutzwiller-projected generalized BCS wavefunctions have been used, finding a gapped Z 2 topological phase throughout; in contrast to this, Ref. 22, supplementing the same ansatz with two Lanczos steps and anisotropic couplings in an enlarged unit cell, finds that around the isotropic point γ = 1, a gapless U (1) Dirac spin liquid (DSL) phase outperforms the gapped Z 2 phase for sufficiently large systems, while for γ 0.25, Valence Bond Crystal (VBC) order dominates. Finally, Ref. 23 analyzes the model using iDMRG, supplemented by exact diagonalization, and finds a U(1) DSL for sufficiently large γ, which at γ 0.1 transitions to a phase with nematic order (i.e., breaking lattice rotation symmetry). In the light of these conflicting results, the nature of the strongly anisotropic limit, and the question whether it is adiabatically connected to the isotropic KHAFM seems wide open.
In this paper, we use an ansatz based on a systematic optimal cooling procedure, applied to the RVB state, to analyze the nature of the breathing KHAFM, focusing on the strong anisotropy limit. Our ansatz, which we term "simplex RVB", clearly outperforms the previous results obtained for the thermodynamic limit, and clearly yields a gapped Z 2 spin liquid rather than a VBC phase. Our ansatz differs from previous approaches in several ways: First, it implements a systematic and optimized cooling procedure -in essence, an optimized imaginary time evo- lution scheme -which can be systematically constructed from any Hamiltonian at hand. Second, it requires only a very small number of parameters with a clear physical interpretation; in our case, we use at most 3 parameters. Third, those parameters have a clear physical interpretation in terms of a variational RVB-type wavefunction: Their role is to create longer-range singlets with suitable amplitude and phase such as to systematically decrease the energy of the variational wavefunction. And lastly, the clear role of the variational parameters in the ansatz facilitates the analysis of its order.
Our analysis reveals a gapped Z 2 topological spin liquid phase for the whole range 0 ≤ γ ≤ 1. In particular, in the strongly anisotropic limit, our results clearly outperform the energies previously obtained in the thermodynamic limit [22] which found a VBC phase, while at the same time they require a significantly smaller number of variational parameters. More specifically, our ansatz with 2 parameters -corresponding to only one optimized trotterized imaginary time evolution step on top of the RVB state -already yields a slightly better energy than the VBC ansatz with 2 Lanczos steps, while it clearly outperforms it with an additional parameter (half a Trotter step). This can be attributed to the fact that our ansatz, unlike Lanczos steps, captures the extensive nature of perturbations and thus correctly reproduces the perturbative expansion in the thermodynamic limit [24].
On a technical level, we use the formalism of Projected Entangled Pair States (PEPS) [25] to implement the simplex RVB ansatz. The idea of the PEPS description is to specify the entanglement structure of the wavefunction as a network of local tensors. The so-called bond dimension determines the efficiency of the PEPS description. The NN RVB state on the kagome lattice can be represented as a PEPS with bond dimension D = 3 [4,26]. While the bond dimension required for p half Trotter steps -corresponding to singlet coverings which contain long-range singlets with range p + 1 -grows as D p = 3 × 2 p−1 , a number p of steps (and thus a singlet span) large enough to yield competitive energies for the anisotropic limit can be reached with computationally accessible bond dimension. A key advantage of this explicit PEPS construction, which is obtained from the RVB PEPS by applying cooling steps, is that it gives us direct access to the relevant entanglement properties for determining the topological nature of the system, and thus allow for a direct and unambiguous identification of the quantum phase of the wavefunction.
The outline of this paper is as follows. In Sec. II we motivate and formally define the simplex RVB ansatz, and give its PEPS construction. In Sec. III, we present our numerical results: First, we discuss the optimal variational energies and corresponding parameters for our ansatz; second, we use PEPS techniques for analyzing the quantum phase of the system as well as the properties of its topological (anyonic) excitations (specifically, anyon masses and order parameters for anyon condensation and deconfinement); and third, we discuss the physical structure of the optimal wavefunction (this is possible due to the clear physical picture behind our variational parameters), as well as possible extensions to potentially further improve the ansatz. Finally, we summarize our results and give an outlook in Sec. IV.

II. SIMPLEX RVB ANSATZ AND PEPS
The key idea behind the simplex RVB ansatz is motivated by the construction used in Ref. 27 to study the isotropic KHAFM. The motivation behind it is to transform a quantum state into the ground state of a given Hamiltonian by algorithmic cooling, that is, by applying local quantum gates (or potentially more broadly local modifications to the wavefunction) which systematically lower the energy where applied. Unlike the construction of the GS in terms of applications of a trotterized imaginary time evolution operator, where all evolution operators are chosen identical and close to the identity, in the simplex ansatz the step size of each evolution operator is optimized variationally such as to minimize the energy. In addition, we start from a well-chosen initial state which already by itself captures essential features of the low-energy physics of the system at hand.
For what follows, it will be convenient to rewrite the Hamiltonian (1) as where P (ijk) is a projector onto the spin-3/2 subspace of 1 2 ⊗ 1 2 ⊗ 1 2 . In order to obtain an approximation of the ground state of H, we can perform imaginary time evolution |ψ = e −βH |φ init for sufficiently large β and a suitable initial state |φ init . If we trotterize e −βH , this yields and accordingly for Q (α). When applied to a suitable initial state, such as the nearest-neighbor RVB state (i.e., a superposition of all nearest-neighbor singlet coverings of the lattice), Eq. (3) has a natural interpretation: First, it is known [28,29] that each NN singlet covering on the kagome lattice contains exactly 25% of "defect triangles", that is, triangles which don't contain a singlet (Fig. 1). Those defect triangles have overlap with the spin-3/2 subspace and thus incur a higher energy than triangles holding a singlet (whose energy is locally optimal). The effect of Q (Q ) is to decrease the weight of defect triangles on right-pointing (left-pointing) triangles. This is achieved by creating longer-range singlets: Since (where R (ijk) rotates the qubits), acting with P (ijk) on a defect triangle produces longer-range singlets, i.e. .
(Note, however, that the linear dependence between pairwise permutations and rotations for 3 qubits -which allows for the form (5) -implies that this way of looking at the longer-range singlet patterns is not unique; a unique pattern within the singlet subspace can be singled out by avoiding crossings.) Unfortunately, representing imaginary time evolution accurately using a Trotter expansion such as Eq. (3) is costly, as it requires a large number of Trotter steps which grows with the system size; in particular, in the context of tensor networks this incurs an exponentially growing bond dimension. One option here is to compress e −βH to a more compact tensor network [30,31]; however, in this process, translational invariance is either lost or obfuscated, and the physical interpretation of the parameters missing. We therefore resort to a different approach: We restrict to a small number of "Trotter" layers in Eq. (3), but we allow for independent parameters α i for each step i and optimize all those parameters such as to minimize the variational energy. This leads to the following ansatz, which we term simplex RVB : where * ∈ { , } is determined by the parity of p, and the α i , i = 1, . . . , p, are the variational parameters. Note that we choose to apply Q leftmost: This way, the ansatz yields the correct behavior in leading order perturbation theory around γ = 0 (we discuss this in detail in Sec. III A); in agreement with this, we observe that this ordering gives better energies in particular in the limit of strong anisotropy.
We now give the PEPS description of the simplex RVB ansatz. We start by reviewing the construction of the NN RVB state [4,26] which is comprised of triangular and on-site tensors. The triangular tensor is defined to be the sum of one configuration with a defect (containing no singlet) and three configurations without defect (containing one singlet each): where δ and ε denote 3-dimensional Kronecker delta and fully antisymmetric tensors, respectively. The on-site tensor ensures that every site is paired with exactly one of its neighbors, The resulting tensor network, obtained by blocking the triangular and on-site tensors, has a three-site unit cell and is given in Fig. 2a. We implement each local action (I − αP ), which is not unitary, as a "controlled" gate on three qubits, controlled by a control qubit (this will be useful for extensions of the ansatz discussed later). The gate acts trivially when the control qubit is |0 , while a projector onto the energetically favorable spin-1/2 subspace of the three qubits is applied if the control qubit is |1 (Fig. 2b). For the time being, we choose the control qubits in a product state |0 + α|1 , leading to a gate Q = (I − αP ), as described previously. For illustration, the tensor network obtained through three applications of Q to the NN RVB, starting with the right-pointing triangles, is shown in Fig. 2c. What is the bond dimension of the simplex RVB PEPS with p layers? We work with the square unit cell shown in Fig. 2a, which contains three kagome spins. With this unit cell and the triangular tensor (7), the NN RVB itself has D = 3. Q (α i ) on right-pointing triangles lie within the unit cell and therefore carry no increase in the bond dimension. Operators Q (α i ) on left-pointing triangles can be implemented with bond dimension 4: E.g., they can be constructed by teleporting the left and bottom neighboring spin to the central site (cf. Fig. 2a), applying Q (α i ), and teleporting them back. D p is therefore multiplied by 4 for every even p, i.e., D p = 3, 12, 12, 48, . . . for p = 1, 2, 3, 4, . . . . However, for p even we can do better: There, the rightmost Q is applied directly to the NN RVB, in which case the state of the teleported spins is already known to the central tensor if the NN RVB index is 0 or 1, allowing to compress the bond dimension for p = 2 to D 2 = 6; thus, we obtain D p = 3 × 2 p−1 = 3, 6, 12, 24, . . . [32].

III. RESULTS
Let us now discuss our results obtained by using simplex RVB states as a variational ansatz for the breathing kagome Heisenberg Hamiltonian (1). The PEPS formalism enables the computation of expectation values of local observables and correlation functions directly in the thermodynamic limit, in contrast to other methods. We use standard numerical methods for infinite PEPS (iPEPS) [33,34], which approximate the boundary by an infinite matrix product state (iMPS) of bond dimension χ (which determines accuracy and computational cost). This allows us to compute the variational energy of an iPEPS with high accuracy and thus to determine the variationally optimal state. In addition, the PEPS approach allows us to utilize the entanglement symmetries of the PEPS and the way in which the iMPS boundary orders relative to those symmetries to study the quantum phase and the topological properties of the optimized wavefunction [35,36]. In our calculations, we choose not to truncate the PEPS tensor before contraction but rather keep the exact simplex ansatz, which avoids truncation errors and gives us direct access to the entanglement symmetries relevant to study the nature of the order in the system; on the other hand, this limits our ansatz to at most p = 3 computationally attainable applications of Q's.

A. Energies
Let us start by giving the results on the optimal variational energy obtained within the simplex RVB ansatz family. For all calculations, we have determined the optimal parameters {α i } through a gradient search using the corner transfer matrix method with a boundary bond dimension χ = 36, and subsequently extracted the energies of the optimized wavefunctions using boundary iMPS (i.e., the fixed point of the transfer operator) with χ = 64 (this only requires truncation along one direction, resulting in a better convergence of the energy in χ). A table with the detailed energies, including a convergence analysis and error bounds, as well as a discussion of a potential extrapolation in p, are given in Appendix A.
In Fig. 3a, we plot the energy density e (i.e., the energy per site) of the optimized simplex RVB wavefunction for the breathing kagome Hamiltonian (1) as a function of the anisotropy γ, for p = 1, 2, 3. For better comparison in the strongly anisotropic limit, we plot in Fig. 3b e(γ) + 0.1353 γ for 0 ≤ γ ≤ 0.2, where the subtracted linear offset corresponds to the behavior in first order perturbation theory, as obtained by extrapolating DMRG calculations [23,37] for the effective first-order model [20] on finite cylinders. Beyond the p = 1, 2, 3 simplex RVB results, we also show the data obtained for the VBC (and the energetically less favorable U(1) DSL) ansatz in Ref. 22 using VMC. We find that already for p = 2, our ansatz gives energies slightly below the VBC ansatz, and for p = 3, it clearly outperforms it. This is particularly remarkable since the p = 2 (p = 3) simplex ansatz has only two (three) parameters, corresponding to effectively one (one and a half) imaginary time evolution steps, while the VMC ansatz of Ref. 22 has 11 parameters, including two Lanczos steps (cf. the discussion in the introduction). In addition, Fig. 3b also shows energies obtained by extrapolating DMRG data for the full model (1) for cylinders with N v = 8, 10, 12 to N v → ∞, which we find to be remarkably close to our p = 3 data in the strong anisotropy regime around γ ≤ 0.04, given that our ansatz only depends on three parameters rather than about 10 8 . Since the extrapolation of the DMRG data is subtle (cf. Appendix B), the gray cone indicates the linear order extracted from the N v = 12 and N v → ∞ DMRG data, which we expect to provide reliable lower and upper bounds to the true slope for the full model.
For better comparison in the strong anisotropy limit,  we expand the energy density for small γ as where the c i depend on p. The values c 1 for the slope at γ = 0 for the different methods are given in Table I (see  Appendix B for details on the extraction). They confirm that for p = 2, our ansatz performs slightly better than the VBC ansatz [22], while for p = 3, it clearly outperforms it. The DMRG data for the nematic spin liquid is the same as the gray cone in Fig. 3b, that is, obtained from the DMRG data of Ref. 23 [37] for N v = 12 cylinders and the N v → ∞ extrapolation, which should give lower and upper bounds to the true value. Finally, we give values for c 1 obtained by extrapolating to p → ∞ in the inverse bond dimension 1/D p ∼ 1/2 p , which we expect to be a reasonable fit in a gapped phase, and which yields a value competitive with the DMRG results. An alternative way to extract c 1 is by using a perturbative expansion. Perturbation theory predicts that for small γ, the energy per site is given by where G denotes the ground state manifold at γ = 0, that is, the subspace of all states with spin 1 2 on the right-pointing triangles. Within our variational family, this corresponds to fixing α 1 = 1 (which explains why we want to have Q leftmost in Eq. (6) if we want to correctly reproduce the perturbative limit), and letting G be the set of simplex RVBs for a given p with fixed α 1 = 1. We thus find that within our ansatz family, we can determine c 1 perturbatively as that is, by minimizing the energy density e on the leftpointing triangles within the simplex RVB family with α 1 = 1. This optimization incurs one less variational parameter and does not require fitting, and can thus be carried out to significantly higher precision; we report the corresponding values in Table I alongside the values obtained from fitting e(γ), which are in excellent agreement.
As an additional check for the quality of the optimal variational state, we have considered the energetics of left-and right-pointing triangles at and around the Heisenberg point. We find that even though the p = 3 ansatz treats the inequivalent triangles differently (in particular, Q acts twice on the the right-and only once on the left-pointing triangles), the energy splitting between the triangles vanishes for the optimal energy wavefunction (Fig. 4a); we observe the same effect also for the optimal wavefunction with p = 2. Alternatively, we can consider the optimal energy density in a symmet- in the vicinity of the Heisenberg point. We obtain a fit e p=2 GS (δ) ≈ −0.4283 + 0.001δ − 0.083δ 2 and e p=3 GS (δ) ≈ −0.4333 − 0.001δ − 0.071δ 2 for the simplex RVB ansatz with p = 2 and p = 3, respectively (Fig. 4b), which essentially show a zero slope at δ = 0 and is thus symmetric around δ = 0 to very good accuracy, as required by symmetry considerations [38].

B. Order, correlations, and quantum phase
The PEPS description of the NN RVB has a Z 2 symmetry on the entanglement degrees of freedom. Such a symmetry has been shown to be essential to explain the topological features of PEPS models, as well as to understand and analyze the breakdown of topological order in such systems [35,36,[39][40][41][42]. This is accomplished by considering the boundary iMPS obtained when contracting the 2D PEPS (i.e., the fixed point of the transfer operator), and analyzing how it orders relative to those symmetries: The specific type of order is directly related to the quantum phase displayed by the bulk wavefunction. From those symmetries, we can construct halfinfinite string operators which on the one hand create anyonic bulk excitations in a given anyon sector a (with a = s, v, sv for spinon, vison, and the composite fermion, respectively), but at the same time form (string) order parameters which detect the ordering of the boundary state. By computing expectations values of these string operators either in one layer (denoted a ) or in both layers (denoted aa † ), we can construct order parameters which probe the condensation and confinement of anyons, and thus the proximity to a topological phase transition; specifically, a non-zero value of the "deconfinement fractions" aa † , as well as "condensate fractions" a ≡ 0, are indicative of the topological phase. At the same time, for vanishing order parameters we can study the rate at which the corresponding expectation value for finite strings with two dual endpoints decays to zero as their separation increases, giving rise to corresponding length scales for condensation (mass gap) and confinement.
We have computed anyonic order parameters for the vison, as well as correlation lengths for all anyons, for the optimized simplex RVB for p = 1, 2, 3 as a function of the anisotropy γ. Since we do not truncate local tensors during optimization, the entanglement symmetries of our tensors remain easily accessible, facilitating the analysis. The corresponding data is shown in Fig. 5. For the anyonic order parameters (Fig. 5a), we find that v = 0 and vv † > 0, which implies that the system is in a Z 2 topologically ordered phase for the given p = 1, 2, 3. The vv † for the different p all show only a small γdependence, with no indication of a phase transition at some intermediate γ building up at larger p. On the other hand, at least for γ close to 1, vv † clearly decreases with p, leaving open the possibility of a critical phase around the Heisenberg point.
Next, let us analyze the correlation lengths, shown in Fig. 5b for p = 3. We find that the dominant correlation length is given by spinon correlations, as known for the NN RVB state [42]. In addition to the different anyon correlations, the figure also shows data for spinspin ( S i .S i+r ) and dimer-dimer ( D i .D i+r ) correlations, computed for the spin and dimer pairs indicated in Fig. 5d. Again, all lengths change smoothly with γ, and while we observe a minor increase of correlations with γ, there is no sign of a phase transition. Note that the similar behavior of spinon and leading trivial (including dimer-dimer) correlations and their relative scale is consistent with previous observations [43] which could be explained as arising from correlations between pairs of spinons [44].
In Fig. 5c, we compare the spinon and dimer correlation lengths ξ p • for the different p = 1, 2, 3. We find a surprising behavior: While the curves for p = 1 and p = 2 show qualitatively similar behavior (with increased correlation length for p = 2), and display a decrease of correlations with growing γ, the p = 3 curve exhibits the opposite behavior. Even more noteworthily, while at the Heisenberg point, correlations keep increasing with p (consistent with a gapless phase), in the small γ regime the correlations decrease again, speaking against a long-range ordered or critical system. The behavior for p = 1, 2 can be qualitatively explained from the way the Q(α) act (cf. the next section where the optimal α are discussed): Q (α 1 ) acts on the strong right triangles. As it decreases the energy of the latter, α 1 will increase with growing anisotropy. At the same time, the Q's create longer-range singlets, which should give rise to longerrange correlations: Correlation functions are obtained from overlaps of singlet configurations, weighted by the number of singlets involved as 2 − /2 , and longer-range singlets allow to connect two points at a given distance with smaller and thus larger weight [45]. The additional increase in correlation length for p = 2 can be explained from the presence of the Q (α 2 ) layer which gives rise to additional longer range singlets before the application of Q (α 1 ). But why does the behavior change for p = 3? As we discuss in more detail in the next section, the role of the topmost Q's is to adjust the energy, as they shift the weight between spin 1 2 and 3 2 right before applying the Hamiltonian; the optimal value of the corresponding α iand thus the amount of correlations they create -is thus governed by immediate energetic considerations (i.e., the overlap with the spin 1 2 space). The lower-lying Q layers, on the other hand, are not directly relevant for the energetics -rather, their job is to set up the underlying wavefunction by creating longer-range singlets in a way where the topmost layers can produce the best possible energies for left-and right-pointing triangles simultaneously. Thus, it is only with the lower layers i ≥ 3 that the Q's primarily serve the purpose of creating the right type of long-range singlets and correlations, rather than just tuning the value of the energy. It therefore seems plausible that the p = 3 behavior of the correlations is closer to the true behavior at large p, and we expect this tendency to continue as we further increase p.
Finally, let us discuss the possibility of a nematically ordered phase, as proposed in Ref. [23] for the strongly anisotropic limit; here, the nematic order was found to break rotational symmetry around the center of the triangles, leading to different Heisenberg energies along inequivalent links. By construction, our ansatz keeps all symmetries of the Hamiltonian and thus cannot explicitly break this symmetry. On the other hand, if nematic order were favored we would expect the system to form a longrange ordered state, that is, an equal weight superposition of all three nematically ordered states. This longrange order should be reflected in a diverging correlation length, and thus, the absence of any such divergence in the dimer-dimer correlations (which in fact rather decrease) speaks against the presence of a nematically ordered phase. To further strengthen this point, we have also considered the dimer-dimer correlations D 0 .D r at a given short distance r, shown in Fig. 5e -if nematic order were favorable, we would expect to see such correlations build up at short distance already at low p. However, Fig. 5e shows no significant increase of D 0 .D r at small γ, and it in fact decreases for p = 3. More importantly, the observed values of D 0 .D r are of the order of 10 −4 even for r = 1, while for the nematic order parameter found in Ref. [23], we would expect them to be on the order of 0.03. Overall, we find that our results show no indications of a nematic phase in the strong anisotropy limit γ 1. As a final test for the Z 2 topological spin liquid nature of the simplex RVB state for p = 3 for all values of γ, we study the deconfinement order parameter for visons and the spinon correlation length as we interpolate from the NN RVB state (which is known to be a gapped Z 2 topological spin liquid) to the optimal simplex RVB wavefunction, by interpolating α(θ) = θ × α γ from θ = 0 (NN RVB) to θ = 1 (optimized simplex RVB), where α γ denotes the optimal parameter values for a given γ. The result is shown in Fig. 6. Again, we find that both quan- tities change smoothly, re-confirming the topological Z 2 spin liquid nature of the optimal wavefunction for the whole range of γ.

C. Structure of optimal wavefunction and possible generalizations
The fact that our simplex RVB ansatz encompasses only very few parameters with a clear interpretation allows us to directly study how the structure of the optimal wavefunction changes as we vary γ and increase p. We recall that our ansatz [Eq. (6)] was of the form where Q • (α i ), • ∈ { , } projects onto the spin-1/2 subspace of the corresponding triangles for α i = 1 and acts trivially for α i = 0 -that is, it lowers the energy of those triangles as α i approaches 1. At the same time, it increases the number of longer-range singlets, as it acts by permuting the singlets; following Eq. (5), one can argue that the largest amount of singlets is permuted at α = 3, though this has to be taken with due care due to the large number of linear dependencies of different long-range singlet patterns, as well as cancellations in the singlet range growth from permutations on adjacent sites. Indeed, the fact that the Q's arise from trotterizing the imaginary time evolution implies that a certain amount of such long-range singlets is required to obtain a good variational wavefunction. Fig. 7 shows the optimal values of {α i } p i=1 for p = 1, 2, 3, as a function of the breathing anisotropy parameter γ. For p = 1, we find that for maximum anisotropy γ = 0, α 1 = 1 -this is expected, as it forces all strong triangles to have spin 1/2 and thus minimum energy, and the state of the weak triangles is irrelevant for γ = 0. As we increase γ, we observe that the value of α 1 decreases, increasing the probability of the weak triangles to have spin 1/2. Remarkably, however, we see that the optimal value of α 1 even at the symmetric point γ = 1 is significantly above 0: This can be understood from the . .  fact that Q in the simplex RVB ansatz does act not only by shifting the weights of defects between inequivalent triangles in the RVB, but at the same time creates energetically favorable longer-range singlets. Note, however, that those singlets are created by decreasing, rather than increasing, α 1 , suggesting that for p = 1, the amount of longer-range singlets at the Heisenberg point is smaller than for γ 1. For p = 2, the behavior of α 1 closely resembles that for p = 1, but with a smaller change (i.e., a larger value α 1 ) towards γ = 1. The correspondingly lower energy gain on the left-pointing triangles is compensated by Q (α 2 ), which lowers the energy of the left-pointing triangles prior to the application of Q (α 1 ), while in addition allowing for the creation of NNN neighbor singlets, thus being energetically favorable.
For p = 3, we make a similar observation: The curves for α 1 and α 2 are again quite close to the p = 2 case. Now, the value of α 2 for small γ has increased -giving the weak left-pointing triangles a lower energy, but also creating longer-range singlets. The energy gain of the weak triangles is now compensated by biasing the system towards strong triangles using α 3 . An interesting point to note is that α 3 > 1, unlike α 1 and α 2 : That is, in the first layer applied to the NN RVB, it is now favorable to flip the sign of the spin-3/2 space or -to the extent the picture of Eq. (5) as creating longer-range singlets is correct -to create a larger fraction of longerrange singlets at the expense of not lowering the energy of the strong triangles as much as possible. Indeed, the latter interpretation is plausible, given that longerrange singlets are overall energetically favorable, and the immediate energetics is taken care of by Q (α 1 ). If we were to follow this reasoning, we would expect further layers to also have α i > 1, i ≥ 3; this suggests that the qualitative change in the behavior of order parameters and correlations occurring at p = 3 will persist for larger values of p.
Given the observation that the lowest layer (i.e., Q • (α p )) significantly biases the NN RVB towards configurations with no defects on one kind of triangles, it seems plausible that a modification of the NN RVB layer in a way which biases it towards configurations with less defects on the suitable triangles should further improve the ansatz. This can be done following the idea of the original simplex RVB paper [27]: We modify the rightpointing triangular tensor (7) of the NN RVB as where a parameter β > 0 (β < 0) effectively shifts the amplitude of defect configurations towards left-pointing (right-pointing) triangles. Importantly, this modification does not lead to an increase in the bond dimension. Subsequently, we can apply the Q's as before, to obtain an enhanced simplex RVB ansatz with one additional parameter. We have tested this ansatz and found that it does not lead to better variational energies, except in the case p = 1. We attribute this to two facts: As discussed, the energetics is predominantly taken care of by the top Q layers (in particular α 1 and α 2 ), rather than the low-lying β; on the other hand, the lower layers mostly serve the purpose to create longer range singlets, whereas β, while changing the weight of the spin-1/2 space on the corresponding triangles, does not give rise to longer range singlets (and in fact reduces the spinon correlation in the system).
As mentioned earlier, we can consider the Q operators as gates which are controlled by a "control qubit" |0 + α i |1 . This enables us to interpret the simplex ansatz as a variational optimization over p control qubits chosen from the manifold of product states (Fig. 2c). In principle, we can further enrich the simplex ansatz (6) by allowing the control qubits to be in a general p-qubit state, i.e., correlating the α i of the different layers. This provides a significantly enlarged simplex manifold, even though it does not allow to increase the range of the singlets. However, we have found that for the computationally feasible values of p, the optimization over the manifold of general p-qubit control states does not lead to an improvement in energy as compared to the optimization over the manifold of p-qubit product states.

IV. SUMMARY AND OUTLOOK
In this paper, we have introduced a simple yet powerful ansatz for the kagome Heisenberg antiferromagnet (KHAFM) with breathing anisotropy, termed simplex RVB. Our ansatz is physically motivated from algorithmic cooling, and effectively consists of p/2 imaginary time evolution layers with optimized step sizes applied to the NN RVB, approaching the true ground state as p → ∞. It yields simple few-parameter families of wavefunctions with a clear physical interpretation in terms of longer range singlets, which are energetically favorable. The ansatz has a simple PEPS representation, which makes it amenable to numerical simulations and an in-depth analysis of its order.
We have analyzed the optimal simplex RVB ansatz for p = 1, 2, 3 for the breathing KHAFM, with a focus on the strong anisotropy limit, and found that already for p = 2 it improves over existing VMC results, while for p = 3 it clearly outperforms them, even though it requires significantly less parameters. We also find that for p = 3, our energies are rather close to extrapolated DMRG energies. It is thus probable that with just a few more layers, the simplex RVB will be fully competitive with DMRG simulations, which is remarkable given the small number of parameters.
We have investigated the nature of the order in the optimized simplex RVB for the breathing Heisenberg model, using a wide range of probes based on the explicit PEPS description and the underlying entanglement symmetries. We find that for the whole parameter regime, our ansatz yields a gapped Z 2 topological spin liquid for the accessible values of p. In the strongly anisotropic regime, we find that the correlations saturate as we increase p, thus exhibiting no signs of long-range order which one would expect e.g. for a nematically ordered phase. On the other hand, at the Heisenberg point, our results show a clear tendency to larger correlation lengths as p increases, which is consistent with a critical DSL phase at the Heisenberg point; both the improvement in energy and the growth of correlations with p points to the relevance of long-range fluctuating singlets for the kagome Heisenberg antiferromagnet.
In order to benchmark a gapped vs. a gapless spin liquid in particular at the Heisenberg point, it would be interesting to compare the simplex RVB ansatz with a variant where one starts from a gapless U(1) spin liquid rather than the gapped NN RVB. One idea which fits well with the PEPS picture is to change the Z 2 invariant PEPS tensors by U(1) invariant ones (which we expect to give a critical wavefunction), e.g. by omitting those Z 2invariant configurations which break U(1). We describe and test such an ansatz in Appendix C. However, while the resulting ansatz indeed yields a gapless spin liquid, we find that the corresponding wavefunction is energetically unfavorable for the Heisenberg model. The reason for this can be found in the general approach of the construction: Since different Z 2 configurations map to different singlet patterns, removing configurations amounts to omitting certain singlet patterns and thus breaks the lattice symmetry, which induces doping with visons and ultimately closes the vison gap. However, as we have observed, the dominating correlation (and thus gap) at the Heisenberg point is given by spinon correlations. Thus, a suitable ansatz would have to drive the system into criticality through doping with spinons. To this end, we would have to resort to a different approach and allow for longer-range singlets e.g. by introducing teleportation bonds in the PEPS [46], which break the Z 2 symmetry. We leave the study of such an ansatz for future work. First, we report in this appendix the energy densities of optimal states within the simplex RVB ansatz for the breathing hamiltonian. Then we discuss the convergence of energies with an increasing bond dimension of environment tensors. In the end, we discuss a possible extrapolation of energies for p → ∞.
We have computed the energy densities for different γ and p for χ ≤ 84. Fig. 8 shows the behavior of the energy densities e χ vs. 1/χ, relative to an arbitrary offset. We see that the fluctuations decrease by roughly an order of magnitude as we increase χ to χ ≥ 60. We therefore choose to estimate the energy and its error by taking the mean and standard deviation of e χ for χ ≥ 60. The resulting energy densities are listed in Table  II. We find that energies tend to converge more quickly in the strongly anisotropic regime as compared to near the Heisenberg point; moreover, their convergence is faster for smaller values of p. Let us note that due to the finite χ, we observe a small splitting in the different bond energies within the left-or right-pointing triangles, of the same magnitude as corresponding flucutations of e χ in Fig. 8. It is instructive to note that the computation of the expectation values of local observables is more expensive in comparison to the calculation of correlation lengths. The calculation of correlation lengths does not require the whole environment, but only the two fixed points of the transfer operator, which enables us to compute them for larger χ's [34,36].
As can also be seen from the data in Table II, the convergence of energies in p is qualitatively different for small and large values of γ ( Fig. 9(a-c)). The energies exhibit negative (positive) curvature in the presence of strong (weak) anisotropy. One could try to argue that the different scaling behavior of energy densities is a manifesta-   The offset e ofs is chosen by averaging eχ for χ ≥ 16. We use eχ for χ ≥ 60 to obtain the mean energies and standard deviations given in Table II, cf. text. tion of the system undergoing a transition from a gapped Z 2 to a gapless SL phase as γ increases, but a meaningful assessment would clearly require further data points in p.
Given that we only have p = 1, 2, 3, and there might well exist even-odd effects at smaller p, extrapolating the energy for p → ∞ is fairly speculative. One attempt would be based on assuming a gapped phase at small γ, which suggests an exponential convergence of the energy in the number of Trotter steps (as states above the gap are exponentially supressed) and thus in p; this amounts to the common extrapolation in the inverse bond dimension 1/D p ∝ 1/2 p . Fig. 9d shows e p=∞ obtained from fitting the odd p with 1/2 p ; the obtained values for small γ are in good agreement with the extrapolated DMRG data (see Table I).
. . In this appendix, we discuss the extraction of the coefficient c 1 of the energy expansion e(γ) = −0.25 + c 1 γ + c 2 γ 2 + . . . in the strongly anisotropic limit. Since c 1 = lim γ→0 ∂e ∂γ , given a sufficient number of points, we can reliably take the numerical derivative of the energy density and examine its behavior in the limit γ → 0. Fig. 10(a-d) shows the derivative of energy density as a function of γ for the simplex RVB. From a linear fit of the data, we can extract c 1 for different values of p. As an alternative approach, we can also directly fit the energy densities quadratically over the interval [0, γ max ], and extract c 1 in the limit γ max → 0. The values of c 1 which we get for the simplex RVB by using either the derivative or quadratic fitting are in very close agreement with each other.
Furthermore, as explained in Sec. III, we can use perturbation theory to extract c 1 for the simplex RVB by fixing α 1 = 1 and minimizing the energy density e ≡ c 1 of left-pointing triangles. The corresponding c 1 are indicated by stars in Fig. 10; they match well with the values we get by fitting the derivative of energies. Fig. 11 shows the energy density of left-pointing triangles for the sim-  .
(f) FIG. 10. (a-d) Numerical derivative of the energy density with respect to γ in the strongly anisotropic regime for the simplex RVB for p = 1, 2, 3, ∞; in (d), the limit p → ∞ has been taken before the derivative. (e) Derivative of the energy density and (f) coefficient c1(γmax) for the DMRG data from Ref. 23 [37]; the limit N → ∞ has been taken before the derivative.
Energy density of left-pointing triangles as the function of parameters α2 and α3 with the constraint α1 = 1 for the p = 3 simpex RVB, used to perturbatively extract c1. Note that the energy landscape possesses no local minima, rendering the optimization stable.
plex RVB with p = 3 as a function of the remaining two parameters; noteworthily, we find that the energy landscape has no local minima and thus allows for a reliable minimization. Let us emphasize here that in the case of simplex RVB, the perturbative approach for evaluating c 1 is more efficient as it eliminates one variable from the variational optimization.
For comparsion, we have extracted the coefficient c 1 for the DMRG data from Ref. 23 [37]. Fig. 10(e) shows the derivative of the energy density for different YCN v -2 cylinders, as well as the derivative of the N v → ∞ energy obtained by linear extrapolation in 1/N v (changing the order of limit and derivative gives almost extactly the same values). However, one needs to be very careful about the extrapolation as the DMRG data suggests a phase transition to the nematic phase in the strongly anisotropic limit, and the phase boundary is sensitive to N v . To extract bounds on c 1 from the DMRG data, we therefore examine the behavior of c 1 (γ max ) which is obtained by linear fitting of ∂e ∂γ over the interval [0, γ max ] for both the N v = 12 and extrapolated N v → ∞ data (Fig. 10f). A further linear extrapolation of c 1 (γ max ), where we restrict to points γ max in the nematic phase, gives estimates of c 1 for the DMRG data.

Appendix C: Simplex ansatz and vison condensed gapless spin liquid states
In RVB-type wavefunctions, U(1) symmetries are typically related to critical behavior. This is most prominent in the square lattice dimer or RVB model, which can be mapped to a height model and thus a critical U(1) field theory [47,48]. Corresponding behavior has also been observed in PEPS, where continuous virtual sym-metries (corresponding to the counting of singlets) have been related to critical behavior, both in the context of RVB-type models and beyond [49][50][51].
In this appendix, we construct a modification of the kagome NN RVB which acquires a U(1) symmetry, and show that it yields a gapless spin liquid state; subsequently, we use this state to build a candidate for a critical simplex RVB and analyze its suitability as a candidate for the breathing KHAFM. For the PEPS description, we begin by considering the blocked local tensor of the NN RVB wavefunction. It is defined on a three-site unit cell containingone right-pointing triangle and is of the form with λ = 1. Here, each term in the sum can be regarded as a tile, marked with a "bond" color (red or black) at each edge. The tiling rule is that only edges with equal color can be connected. If we subsequently ignore the black edges, this results in a pattern where pairs of neighboring vertices on the kagome lattice are connected by red lines: If we replace those by singlets, this precisely yields a NN singlet covering, and it is straightforward to check that all NN coverings areobtained this way. (At a technical level, the bond dimension is D = 3, where black legs correspond to a |2 "no singlet" state, while red legs are in the states |0 , |1 , and perfectly correlated with the corresponding physical spin. Pairs of red vertices within a unit cell are placed in a singlet, and contraction involves an iσ y in the {|0 , |1 } subspace, see Refs. 4 and 26).
The red legs clearly obey a Z 2 Gauss law with a −1 charge per unit cell. We can now introduce a generalization of the NN RVB by changing λ, which corresponds to changing the relative weight of different singlet coverings. In particular, if we choose λ ≡ 0, we are left only with terms which all have exactly one red leg: The tensor A(λ ≡ 0) has acquired a (staggered) U(1) Gauss law, in analogy to the square lattice RVB state, suggestive of a critical behavior at γ = 0.
We start by analyzing the behavior of the "bare" RVB as we change λ. Fig. 12a shows the correlation lengths for trivial and anyon-anyon correlations for the different anyon sectors. We find that as we approach λ → 0, the vison and trivial correlations diverge; this is consistent with the interpretation that changing λ, which introduces a disbalance between different singlet configurations, dopes the systems with visons which ultimately makes the vison gap close and gives rise to criticality due to vison condensation. The critical nature at λ = 0 is confirmed by analyzing the dependence of the correlation length at λ = 0 on χ, which exhibits a rapid growth. Note that we observe criticality only at λ = 0, which is quite different from the model in Ref. [52] where an extended critical region around the U(1) point is reported. We now use the U(1) RVB (i.e., at λ = 0) to construct a simplex RVB, that is, we follow Eq. (6), but with the λ = 0 RVB as a starting state. We have computed the optimal energies for p = 3, shown in Fig. 12b. However, we find that the energies are above the energies for the simplex RVB ansatz built on top of the NN RVB state, in particular in the vicinity of the Heisenberg point. This is consistent with the fact that we expect a closing of the spinon gap to be the driving mechanism behind a potential critical spin liquid at the Heisenberg point, while our ansatz at λ = 0 exhibits criticality to due a closing vison gap. For the same reason, the application of Q's tends to increase the spinon correlations at the cost of decreasing vison correlations, thus effectively driving the ansatz away from criticality (at least initially); this can also be understood from the fact that the Q's permute singlets and thus can restore the singlet patterns which are missing in the RVB at λ = 0.