Spin-$\frac{1}{2}$ kagome Heisenberg antiferromagnet with strong breathing anisotropy

We study the zero-temperature phase diagram of the spin-$\frac{1}{2}$ Heisenberg model with breathing anisotropy (i.e., with different coupling strength on the upward and downward triangles) on the kagome lattice. Our study relies on large scale tensor network simulations based on infinite projected entangled-pair state and infinite projected entangled-simplex state methods adapted to the kagome lattice. Our energy analysis suggests that the U(1) algebraic quantum spin-liquid (QSL) ground-state of the isotropic Heisenberg model is stable up to very large breathing anisotropy until it breaks down to a critical lattice-nematic phase that breaks rotational symmetry in real space through a first-order quantum phase transition. Our results also provide further insight into the recent experiment on vanadium oxyfluoride compounds which has been shown to be relevant platforms for realizing QSL in the presence of breathing anisotropy.


Introduction
Quantum spin-liquids [1] are exotic phases of matter with highly entangled ground-states and fractionalized excitations which fall beyond the Ginzburg-Landau paradigm [2][3][4]. Due to lack of local magnetic ordering at zero temperature, QSL are very hard to come by experimentally [1]. However during the past years, there has been tremendous progress at the theoretical level towards a better understanding of these phases, ranging from the QSL ground-states of topological quantum codes [5][6][7][8][9], which are platforms for fault-tolerant quantum computation [5,6], to highly frustrated quantum spin systems such as Heisenberg antiferromagnets on lattices with triangular geometries [1,[10][11][12][13].
The AFHK is also of experimental relevance. Recent experiments suggested candidate materials, such as herbertsmithite [33], volborthite [34,35], or vesignieite [36] for realizing AFHK. More importantly for our present purpose, it has been suggested that kagome compounds with breathing anisotropy, i.e., different Heisenberg exchange coupling on the upward (J ) and downward triangles (J ) can still host QSL groundstates even at relatively large breathing anisotropy [37]. Signatures of a gapless spinliquid have recently been observed experimentally in a vanadium oxyfluoride compound, [37][38][39], at J /J ≈ 0.55. This finding is of particular importance since naturally kagome compound tend to be anisotropic in nature due to impurities or effective perturbations which influence the strength of coupling on different kagome triangles. The spin- 1 2 breathing-kagome Heisenberg model (BKH) of Eq.(1) has first been studied in the nineties [14,40]. An effective model has been derived in terms of two local degrees of freedom per strong triangle, and a mean-field decoupling of these degrees of freedom has led to the conclusion that a family of short-range resonating valence bond (RVB) states could provide a good variational basis [14,41]. Following up on the recent experimental interest in that system, several theoretical studies [42][43][44] aiming at providing further insight into the experimental results have revisited that model. While an earlier VMC study predicted a stabilized Z 2 gapped spin-liquid ground-state for the BKH model for J /J ≈ 0.7 [43], other studies based on DMRG on semi-infinite cylinders [44] and more refined VMC calculations [42] predicted the persistence of the gapless U(1) spin-liquid up to large breathing anisotropy. The DMRG results further provide evidence that the BKH model undergoes a phase transition to a nematic phase at J /J ≈ 0.2 or below [44] while the VMC results are in favor of a valence-bond-crystal (VBC) which is stabilized for J /J 0.25 [42]. Note that, in this parameter range, a simplex Z 2 RVB state represented as a simple TN was found to have a slightly higher energy [42].
Motivated by the importance of recent experiment on vanadium compounds [37][38][39] and the quest for experimental realization of QSL phases, we address in this paper the spin-1 2 breathing-kagome Heisenberg problem and study, in detail, the full phase diagram and ground-state properties of the BKH model in the thermodynamic limit, by resorting to large-scale tensor network calculations based on infinite projected entangled-pair state (iPEPS) [45][46][47][48] and infinite projected entangled-simplex state (PESS) [28,32] methods on the infinite 2D kagome lattice. In particular, we focus more on the large breathing anisotropy limit and perform accurate energy analysis and finite-size entanglement scaling of energies to try and reveal the true ground-state of the system out of the energetically competing VBC, Z 2 QSL and nematic phases. Our results suggest that the U(1) QSL phase of the isotropic kagome Heisenberg antiferromagnet is stable up to very large breathing anisotropy J /J ≈ 0.05 and that, for larger anisotropy, it undergoes a firstorder quantum phase transition (QPT) to a critical lattice-nematic phase. We capture the lattice-nematic ordering by accurate analysis of the energy density on every bond of the up and down triangles of the kagome lattice in translationally invariant unit-cells with different sizes and further reveal the critical nature of the lattice-nematic phase showing in particular power-law spin-spin correlations along the emerging chains.
The paper is organized as follows. In Sec. 2 we introduce the BKH model on the kagome lattice and briefly discuss the details of the iPEPS and PESS machinery we used for evaluating the ground-state of the system. Next, in Sec. 3 we elaborate on the the U(1) gapless QSL phase of BKH model at the isotropic point with no breathing anisotropy. We further study the large breathing anisotropic limit of the BKH model and lattice-nematic ground-state of the system in Sec. 4. In Sec. 5 we investigate the quantum phase transition and full phase diagram of the BKH model. Finally Sec. 6 is devoted to a discussion and to a conclusion.

Model and Method
The spin-1 2 breathing-kagome Heisenberg antiferromagnetic model [14,40] is defined by where the first (second) sum runs over edges of the upward, , (downward, ) triangles of the kagome lattice (see Fig. 1). Here J and J are the antiferromagnetic (AF) Heisenberg exchange couplings, respectively on the up and down triangles and S i is the spin operator at lattice site i. As discussed above, we are interested in analyzing the ground-state Figure 2: (Color online) (a) Coarse-graining the kagome lattice to a square network of block sites with three spin-1 2 in each block and a total local Hilbert space of d = 2 3 = 8. The iPEPS simulations were performed on 2 × 2 and 4 × 4 unit-cells of block-sites with 12 and 48 spins, respectively. (b) Coarse-graining of a PESS with three physical tensors and two simplices into a block-site. (c) Coarse-graining of a PESS with six physical tensors and four simplices into block-sites.
properties of Hamiltonian (1) for the full ranges of the J /J coupling. Here without loss of generalities, we set J = 1.
In order to study the ground-states of the BKH model, we used two different tensor network methods based on infinite projected entangled-pair state [45][46][47][48] and infinite projected entangled-simplex state [28,32]. The iPEPS algorithm was applied to the BKH problem by coarse-graining the kagome lattice to a square TN of block-sites (each block with three spin- 1 2 ) of rank-5 tensors with physical bond dimension d = 2 3 = 8 and virtual bond dimension D [13,[49][50][51]. The coarse-grained lattice is illustrated in Fig. 2-(a). In order to approximate the ground-state of the BKH model, we used both simple-update (SU) [26,50,52] and full-update (FU) [13,[53][54][55] based on imaginary-time evolution (ITE) [47,56] accompanied by corner transfer matrix renormalization group (CTMRG) [48,52,55] for the contraction of the infinite 2D lattice, applied to 2 × 2 and 4 × 4 unit-cells of block-sites with 12 and 48 spins, respectively. The PESS, however, was performed using 3PESS machinery [32] with rank-3 simplex tensors with no physical leg, which captures the three-partite entanglement inside kagome triangles, and rank-3 physical tensors representing local spins. The contraction of the 2D PESS tensor network was performed by first coarse-graining two simplices and three physical tensors into a block-site with local Hilbert space of d = 2 3 = 8 and then using the CTMRG. We used two different unit-cells for the PESS simulations, one with three physical tensors and two simplices ( Fig. 2-(b)) and another one with six physical tensors and four simplices ( Fig. 2-(c)). While both unit-cells are suitable for capturing U(1) QSL, Z 2 RVB state and lattice-nematic phases, the VBC ordering is only commensurate with the latter unit-cell [42]. We further used simple update based on ITE and high-order singular value decomposition (HOSVD) [32] to optimize the PESS tensors.
The ITE iterations for both iPEPS and PESS methods were performed starting from δτ = 10 −1 down to 10 −5 with at least 4000 iterations for each δτ . We further initialized the simulations with different initial tensors ranging from totally random states to a D = 3 RVB state for different regimes of the problem. Moreover, for each bond dimension D, we recycled both ground-state tensors and the environment from the previous lower bond. In order to check the convergence of the algorithms, we compared the energy difference of two consecutive iterations with a tolerance threshold of 10 −14 for each ITE δτ . In the end, the expectation values of local operators and variational energies were calculated by using the converged ground-state tensors and full environment tensors obtained from CTMRG Refs. [11,25,44].
algorithms [52,54]. Let us further note that within our available resources and time, we managed to push the iPEPS simulations up to D = 6 for the full-update and D = 11 for the simple-update. The PESS simulations were also performed up to D = 13 with the full calculation of the environment, i.e., with CTMRG and up to D = 18 for the simple environment, i.e., with λ matrices of the HOSVD optimization. However, since for the simple environment the energies are no longer variational we do not show them here. Nevertheless, the results were consistent with full environment calculations. Moreover, we fixed the boundary dimension to χ = 80 for D ≤ 8 and χ = 64 for D > 8 in both the iPEPS and PESS simulations. Whenever possible, we further checked the convergence of the simulations for χ > 64 for various D values, and we have noticed no significant change in energies for χ > 64 in our results.
In order to avoid confusion between different results obtained from different methods, we fix the following style in the figures throughout the paper: The iPEPS data obtained with simple-(full-) update are demonstrated with circles (squares) and the PESS data are depicted with diamonds. We have further used stars for illustrating non-variational PESS energies which are obtained from λ matrices of the HOSVD optimization.
In the following, we analyze the BKH Hamiltonian (1) for all ranges of AF Heisenberg exchange coupling and map out the full phase diagram of the system. More precisely, we calculate the ground-state of the BKH model for 0 ≤ J /J ≤ 1 and detect possible phase transitions by following the ground-state energy per site, ε 0 , and its derivative, as well as two-point spin correlators. However, before we start analyzing the phase diagram of the BKH model, let us first consider the extreme limits of Hamiltonian (1).

Small Breathing Anisotropy
In the limit of the BKH model where J /J = 1, i.e., the isotropic point, it has already been shown that tensor network simulations favor a gapless U(1) QSL with algebraic decay of correlations. We refer the interested reader to Refs. [27][28][29]32] for detailed discussions on the isotropic spin-1 2 kagome Heisenberg antiferromagnet. Nevertheless, in order to check the accuracy and efficiency of our algorithms, we applied both iPEPS and PESS methods to the BKH model at the isotropic point and recovered the results of Ref. [28,32]. Fig. 3 illustrates finite-entanglement scaling of the iPEPS and PESS energies for both simple and full calculation of the environment with respect to bond dimension D. We find that ε 0 converges algebraically with D with power-law form ε 0 (D) = ε 0 + aD −β , as shown in the inset of the figure for both iPEPS and PESS energies. The powerlaw convergence of energy with respect to bond dimension has already been suggested as a good numerical signature of critical (gapless) systems, both in one [56] and two [27,28] dimension. Our energies obtained with both iPEPS and particularly with PESS shows the algebraic convergence, consistent with a gapless ground state. Besides, we have contrasted the power-law decay of energies, ε 0 (D) = ε 0 + aD −β , against exponential, ε 0 (D) = ε 0 +exp(−D) and logarithmic decay, ε 0 (D) = ε 0 +log(D) (not shown in the paper) and found that the converges of our energies for both PESS and iPEPS simulations are best fitted with a power-law decay signaling a gapless underlying state.
We have further analyzed the energies on the upward and downward triangles of the kagome lattice as well as, correlations S α i .S α j (α = x, y, z) on every link of each triangle in the unit-cell and and found the same energies and correlations on both triangles with a uniform ground-state which respects all lattice symmetries and SU(2) symmetry of the Hamiltonian (1), thus indicating the QSL nature of the ground-state at the isotropic point.
Our best variational energy at the isotropic point is ε 0 = −0.436979, obtained from PESS with D = 13 which is lower than that of Ref. [32] for D = 13 9-PESS. Let us further note that our iPEPS energy at this point is ε 0 = −0.433374 for D = 11, which is slightly higher compared to that of the PESS due to the lower maximum achievable bond dimension. Besides, it is already known that the simplex structure of the PESS tensors can capture the three-partite entanglement inside a frustrated kagome triangle better and, hence, is able to yield lower energies compared to iPEPS which only captures bipartite entanglement on bonds of the lattice [28,32]. Let us further note that characterizing the U(1) or Z 2 nature of the state is a numerically nontrivial task. The gapless spin liquid is expected to have long-ranged entanglement and correlation functions and the U(1) state has no well-characterized topological order. We calculated the reduced density matrix of the ground state obtained from our PESS simulation on bipartitions of a semi-infinite cylinder and calculated the entanglement entropy [57]. We found no topological contribution to the entanglement entropy which is in agreement with that of a U(1) gapless spin liquid. Moreover, the momentum-dependent excitation spectrum which was extracted from the DMRG transfer matrix of Ref. [29,44], exhibits Dirac cones that match those of a π-flux free-fermion model (the parton meanfield ansatz of a U(1) Dirac spin liquid). Although we are not able to calculate the same information within the framework of our TN simulations, we provided evidences through the paper that our results are in agreement with these DMRG calculations of Ref. [44] and therefore we believe that the gapless state we capture at the isotropic point should be a U(1) spin liquid.
Immediately away from the isotropic point, i.e., for J /J < 1, the difference in the J and J couplings in Hamiltonian (1) breaks the lattice inversion center (i.e. the 180 o rotational symmetry). However, as long as the breathing anisotropy is not too large, we expect the QSL to remain stable. Our analysis indeed confirms that introducing small breathing anisotropy does not destroy the QSL ground state of the spin-1 2 kagome Heisenberg antiferromagnet and the uniformity and SU(2) invariance of the ground-state are preserved at the level of each individual triangles. Nevertheless, upward triangles will have lower energies due to larger J couplings compared to J . Fig. 4 shows the Table 1: Variational energies obtained from iPEPS and PESS compared with U(1), Z 2 QSL and VBC energies of Ref. [42] and DMRG data of Ref. [44,58] for YC8-2 cylinder geometry. Z 2 QSL * is an improved RVB ansatz [59]. The lowest possible energies for all couplings belong first to the DMRG results and next to our PESS (D = 13) simulations. ground-state energy of the system for both iPEPS and PESS at J /J = 0.3. The power law scaling of energies with respect to D is similar to the isotropic point, J /J = 1, which suggests that both points belong to the same phase. In fact, in future sections, and especially in Fig. 8, we will show that the U(1) QSL ground-state of the BKH model at the isotropic point persists to very large breathing anisotropies. We postpone further discussions regarding the stability and persistence of the U(1) QSL ground-state of the system in the presence of breathing anisotropy to Sec. 5.

Large Breathing Anisotropy
In this section, we elaborate on the less studied limit J /J 1, i.e. the large breathing anisotropic regime. In the extreme case where the couplings on the down triangles are zero, J = 0, the system is composed of decoupled upward triangles with AF interactions, with a highly-degenerate ground-state and a ground-state energy ε 0 = −0.25.
By switching on and gradually increasing the J couplings on the downward triangles, forming an ordered ground-state becomes a highly non-trivial task. An early study based on the short-range RVB basis [41] suggests that a gapped Z 2 spin-liquid ground-state may emerge in the presence of large breathing anisotropy. In another study based on Gutzwiller projected wave functions, the analysis of the energy shows that the U(1) spinliquid phase of the BKH model undergoes a dimer instability at large breathing anisotropy and a VBC phase is stabilized for J /J < 0.25 [42]. By contrast, recent DMRG study on semi-infinite cylinders suggest a lattice-nematic ordering for the ground-state of the BKH model at large breathing anisotropies [44] which emerges around J /J ≈ 0.14.
In order to further investigate the nature of the ordering of the system at large breathing anisotropies, we performed systematic energy analysis for the whole range of couplings. In particular, we performed accurate entanglement-scaling of energies versus inverse bond dimension, 1/D up to D = 13 for several points at very large anisotropies. Fig. 5 depicts the scaling of the ground-state energy versus inverse bond dimension, 1/D, for J /J = 0.01. Our energies obtained with both full-update (FU) and simple-update (SU) of the iPEPS as well as PESS results with full calculation of the environment are reported in Table 1 and compared to the results of Ref. [42,59] for U(1), Z 2 QSL and VBC as well as the DMRG results of Ref. [44,58]. Our best variational energy is ε 0 = −0.251363 which was obtained from PESS (D = 13). Our energies are also in agreement with DMRG results of Ref. [44,58] which is slightly lower than our iPEPS and PESS results.
Next, in order to identify the nature of the instability in the large breathing anisotropy limit, we analyzed the energy density and nearest-neighbor spin-spin correlations on every bond of the upward and downward triangles of the kagome unit-cells (see Fig. 2) with different structures. We observed that energies and correlations on different links of the up and down triangles of the kagome lattice undergo a dimer instability with broken C 3 rotational symmetry at the level of individual triangles which are consistent with the lattice-nematic pattern. Fig. 6 illustrates the lattice-nematic pattern with strong and weak correlation on the bonds of the kagome lattice that we captured within the framework of our TN simulations in the large breathing anisotropy limit at J /J = 0.01 for all bond dimensions D and unit-cell structures. Let us further point out that depending on the choice of the initial states, size of the unit-cells and the virtual bond dimension, we observed a spontaneous breaking of mirror symmetry on the up triangles in our TN simulations. The broken mirror symmetry was detected from the difference in the correlations on the two green links of up triangles, 1 , 2 , shown in Fig. 6. We observed that this difference in the correlations of the 1 , 2 links decreases upon increasing the bond dimension, and that the mirror symmetry is restored in the thermodynamic limit. The same behaviour has already been observed in the DMRG simulations of Ref. [44,58].
The nematic state breaks the C 3 local rotational symmetry of the triangles while preserving the translational invariance in every direction of the lattice. This is in contrast with the VBC phase which breaks both translational and rotational symmetries of the system [42]. The lattice-nematic state on the kagome lattice is, in fact, three-fold degenerate. Repeating the simulations with different initial states, we found the two other degenerate nematic states with the same magnitude of correlations on bonds but a different pattern. This implies that in the regime of very large breathing anisotropy, the system undergoes a dimensional reduction with three degenerate ground states that consist of almost decoupled chains. Let us further stress that our results for the large breathing anisotropies are Table 2: The c 1 , c 2 coefficients of the Taylor expansion for the BKH model at large breathing limit obtained with PESS (D = 13) compared with those of the DMRG data for the effective model in Ref. [44,58] (after extrapolation to infinite cylinders) as well as with the Z 2 , U(1) QSLs and VBC states of Ref. [42] and Z 2 QSL * of Ref. [59].
Wave Function [42] −0.1190 −0.079 Z 2 QSL [42] −0.1245 0 VBC [42] −0.1255 −0.055 Z 2 QSL * [59] −0.1323 −0.0628 Effective Model [44] −0.1353 0 in agreement with the recent DMRG results [44]. For completeness of our analysis and further to see how the lattice-nematic state competes with other states from previous studies, we consider the Taylor expansion of the BKH energy in the large breathing anisotropy limit up to second order in perturbation theory The first constant term in the above equation is the energy of decoupled upward triangles. The coefficients c 1 , c 2 can further be obtained by quadratic fits of the energy curve. Table 2 provides the c 1 , c 2 expansion coefficients of our lattice-nematic states obtained with PESS (D = 13) compared with those of the effective model known as trimerized kagome model, which corresponds to a frustrated spin-orbital model on the triangular lattice [14,44]. Note that, since our energy curve shows a discontinuity in its slope, we have performed two separate fits on each side of the discontinuity. The coefficients for Z 2 , U(1) QSLs and VBC state of Ref. [42] and Z 2 QSL * of Ref. [59] as well as those of the nematic state obtained with DMRG in Ref. [44] are also provided in the table. One can clearly see that the c 1 coefficient of the nematic state obtained both in our simulations and previous DMRG results is larger (in absolute value) than those of the Z 2 , U(1) QSLs and VBC states suggesting a stabilized nematic state as the true ground state of the BKH model in the large breathing anisotropy limit.
We have further calculated the long-range spin-spin correlation C(r) = S (x,y) .S (x+r,y) − S (x,y) . S (x+r,y) in the large breathing anisotropy limit compared with several points in the small-breathing regime. Fig. 7-(a) shows the log-log plot of C(r) obtained with PESS for different bond dimensions at the isotropic point, J /J = 1. The C(r) at small distances shows a power-law behaviour with a tail which deviates from the power-law fit (dashed line) at large distances. However, by increasing the bond dimension, D, the spin-spin correlation tends to approach the power-law fit indicating a gapless ground-state for the isotropic point, at the thermodynamic limit. Fig. 7-(b) further demonstrates the C(r) obtained with PESS (D = 11) at J /J = 0.01 in the lattice-nematic phase and J /J = 0.3 compared with the isotropic point in the spin-liquid phase. C(r) in the QSL phase behaves similarly in different directions of the lattice, as expected from a QSL phase with no broken symmetry. Most importantly, C(r) for J /J = 0.3 and 1 decay similarly especially at small distances, suggesting that they belong to the same phase. This is another strong signature that the QSL phase persists in the large breathing anisotropy limit. However, C(r) in the lattice-nematic phase is different along the strong chain, x-direction, and perpendicular to the chain, y-direction, as shown in Fig. 7-(b). We have therefore approximated the power-law decay of C(r) at J /J = 0.3 and 1 and J /J = 0.01 along the strong chain with C(r) ∼ C 0 r −α , shown as a dashed blue line in the Fig. 7-(b). C(r) decays algebraically with C 0 = 0.18, α = 2.2 indicating that the wave functions at J /J = 0.3 and 1 and at J /J = 0.01 along the strong chain are critical. The spin-spin correlation perpendicular to the strong chain at J /J = 0.01 shows a different power-law decay with C 0 = 0.18, α = 3.5, which is the expected behavior of the lattice-nematic state with a different pattern in the x-and y-direction. Note that our findings for the large breathing anisotropy limit are in agreement with the Lieb-Schultz-Mattis theorem [60,61] which states that, for systems with half-odd integer spins in the unit-cell, there cannot exist a gapped spin-liquid with a unique ground state. Therefore, the critical nature found in both phases in the phase diagram is consistent with the theorem.
Let us note that our TN ansatz suffers from a small spurious magnetic ordering of the order ∼ 0.003 which is believed to be an artifact of the PEPS methods with finite bond dimension. The effects of such a spurious magnetization is then seen as zigzag oscillations in the S (x,y) .S (x+r,y) . We, therefore, subtract the local magnetic contribution from the spin-spin correlation to correctly capture the power-law decay.

Quantum Phase Transition
In previous sections, we have identified and characterized two different phases at the extreme regimes of the BKH Hamiltonian (1), the U(1) spin-liquid at the isotropic point J /J = 1 with zero anisotropy and the lattice-nematic phase at large breathing anisotropy limit, J /J 1. It is, therefore, reasonable to expect a quantum phase transition between the two extreme phases. In order to study the QPT, we analyzed the whole regime In an attempt to draw conclusions about the infinite D limit, we have kept track of the evolution of the coefficients c 1 and c 2 with 1/D. First of all, we have estimated the error bars on these coefficients. Fig. 9-(a-d) shows the fitting of PESS energies (D = 13) using the c 1 and c 2 coefficients of Table 2. The error bars on the energies, estimated from the standard deviation from the mean energy values, are very small, of the order 10 −5 − 10 −6 for different couplings, but they still lead to a non negligible error threshold of the order 10 −3 on the c 1 and c 2 coefficients. The resulting estimates for c 1 are plotted as a function of 1/D in Fig. 10, together with the DMRG results of Ref. [44,58] for the effective first-order model in the left panel, and with those of the RVB wave function of Ref. [59] in the right panel. Several remarks are in order. First of all, the values of c 1 for the nematic PESS and for DMRG are very close to each other, a nice confirmation that these approaches are indeed describing the same phase. However, this plot also shows the difficulty in extrapolating our results as a function of 1/D. The naive linear extrapolation suggested by the data if one forgets about the error bars can be excluded because it would lead to a value of c 1 much too small as compared to DMRG, and taking the error bars into account leads to a very broad distribution. Besides, there is no theoretical prediction for the evolution of the results as a function of 1/D on which to rely to go beyond the naive linear extrapolation. The comparison with DMRG suggests that the coefficients c 1 must level off at larger D, but we are unfortunately unable to check this. Finally, we note that the coefficient c 1 of the p-RVB wave-function is not far but above the PESS values  for both the nematic and U(1) phases, in agreement with our finding that, at least for finite bond dimension D, this phase is not stabilized. In order to further check the stability of the U(1) QSL on the right side of the QPT point, we performed a careful entanglement scaling of energies for several couplings explicitly for regions previously suggested to host a phase transition, and we investigated the energy densities and bond correlations on all bonds of the upward and downward kagome triangle. Fig. 11 shows the correlations on all links of the up and down kagome triangles for 0 J /J 1. One can clearly observe that, in the range 0.1 J /J 0.5, the correlations on the links of upwards and on the links of downward triangles are the same, indicating that the C 3 rotational symmetry is preserved at the level of each individual triangle. Besides, we performed entanglement scaling of energies for all of the points in Fig. 11 and observed algebraic decay of correlation for all points. Our findings therefore suggest that the algebraic U(1) spin-liquid phase of the spin-1 2 kagome Heisenberg antiferromagnet is stable up to very large breathing anisotropies. Our results are in agreement with the recent experiments on vanadium oxyfluoride compounds [37][38][39] which detected signatures of a gapless U(1) QSL at breathing ratio J /J ≈ 0.55 [37].

Discussion and outlook
First introduced as a toy model to study the spin-1/2 kagome antiferromagnet [14], the breathing spin-1/2 kagome antiferromagnet has recently attracted attention on its own due to its experimental relevance for some vanadium compounds [37][38][39]. In the present paper, we have performed large scale tensor network calculations of that model based on projected entangled-pair state and projected entangled-simplex state methods. The picture emerging from these calculations is consistent with the DMRG results of Ref. [44]: The system seems to be a U (1) liquid from the isotropic limit J /J = 1 down to very small values of J /J , and to undergo a first-order transition into a critical, lattice nematic phase that breaks rotational symmetry in real space. Our simulations for the largest available values of the bond dimension D locate quite convincingly the transition at J /J 0.05. This should be contrasted to the conclusions of Ref. [44], where the actual critical value of J /J (if any) could not be pinned down because of the still strong dependence of the results on the circumference of the cylinder. However it has not been possible to extrapolate the results to infinite D, and strictly speaking one cannot exclude that the U(1) spin liquid remains stable below J /J 0.05 on the basis of the present results.
In view of the competing phases that have been proposed over the years for the strong breathing limit, let us conclude with a critical review of the situation. The early mean-field theory of the effective model [14] had identified the short-range RVB basis as a promising variational basis, a conclusion partly confirmed in Ref. [41]. Interestingly, recent results based on VMC [42] favour a VBC configuration that is consistent with a particular member of the short-range RVB basis in which all pairs of triangles coupled by a singlet bond are parallel to a given direction. Such a state breaks the rotational symmetry in real space, but it is different from the nematic state reported here because it breaks an additional mirror symmetry.
The results of the present paper and of Ref. [44] on one hand, and the results of a very recent investigation based on an improved simplex RVB PEPS ansatz on the other hand [59], all suggest that energy can be further gained with respect to this VBC state by additional quantum fluctuations, but in different ways.
In the lattice-nematic state, these fluctuations take place along the chains of the phase with broken rotational symmetry. The additional mirror (or equivalently unit lattice translation) that is broken in the VBC state is restored, but the rotational symmetry is not. The ground state consists of very weakly coupled chains. On each chain, the wavefunction is typical of a spin-1/2 critical chain, with correlations that require to go beyond the short-range RVB basis and to include longer-range singlets. So the energy gain w.r.t. the VBC has to do with the development of long-range correlations.
By contrast, in the RVB phase, the full symmetry is restored, the energy gain comes from resonances between short-range RVB configurations, and the spectrum is gapped.
Quite interestingly, the energies of these two states are very close, and even if that of the RVB spin liquid is higher for D = 12 1 , it is not possible to make a reliable prediction for what happens in the infinite D limit. So, if it seems very likely that there is a phase transition at strong breathing anisotropy, the actual nature of the ground state in that limit, a critical lattice-nematic state or a gapped RVB state, should still be considered as an open question. This is a motivation to further improve tensor-network algorithms, and work is in progress along these lines.
In any case, the conclusion that the transition takes place at a very small value of J /J , of the order of 0.05, is in agreement with the recent experiments on vanadium compounds [37][38][39] in which evidence for a gapless U(1) ground-state was observed on a kagome lattice with breathing anisotropy at J /J ≈ 0.55.
Finally, we would like to note that the TN algorithm that we have developed for the kagome lattice can also be used for further investigations of the BKH model in the presence of a magnetic field to capture possible magnetization plateaus and their underlying exotic phases. This is left for future investigation.