Critical colored-RVB states in the frustrated quantum Heisenberg model on the square lattice

TL;DR

This paper introduces colored-RVB states within PEPS to model the quantum disordered phase of the frustrated J1-J2 Heisenberg model, revealing critical correlations near a potential RK point.

Contribution

It constructs a family of SU(2)-symmetric PEPS for the J1-J2 model, showing critical dimer correlations and short spin correlations near an RK point.

Findings

Critical algebraic dimer-dimer correlations at J2/J1~0.55

Short spin-spin correlation length in the same regime

Evidence for a spin liquid or weak Valence Bond Crystal phase

Abstract

We consider a family of SU(2)-symmetric Projected Entangled Paired States (PEPS) on the square lattice, defining colored-Resonating Valence Bond (RVB) states, to describe the quantum disordered phase of the $J_1-J_2$ frustrated Heisenberg model.For $J_2/J_1\sim 0.55$ we show the emergence of critical (algebraic) dimer-dimer correlations -- typical of Rokhsar-Kivelson (RK) points of quantum dimer models on bipartite lattices -- while, simultaneously, the spin-spin correlation length remains short. Our findings are consistent with a spin liquid or a weak Valence Bond Crystal in the neighborhood of an RK point.

Figures

Equations

$$H=\sum _{i,j}$$

$$a=\sum _{\alpha =1}^{\mathcal{D}}{c}_{\alpha }{t}_{\alpha }.$$

$$V=\stackrel{\mathcal{N}\mathrm{times}}{\overbrace{\frac{1}{2}\oplus \cdots \oplus \frac{1}{2}}}\oplus 0,$$

$${\mathcal{G}}_{\beta }\equiv \frac{\partial E}{\partial c_{\beta }}\simeq \frac{E\left[\{c_{\alpha }{\}}_{{}_{\beta }}\right]-E\left[\{c_{\alpha }\}\right]}{\delta }$$

$$C_{\rm d}(r)=\big{<}D_{\bf x}D_{{\bf x}+r{\bf e}_{x}}\big{>}-\big{<}D_{\bf x}\big{>}\big{<}D_{{\bf x}+r{\bf e}_{x}}\big{>}\,,$$

$${\xi }_d(\chi )\sim a_D\frac{\chi }{D^2}+b_D,$$

$$C_{\mathrm{d}}(r)\sim r^{-{\alpha }_{\mathrm{d}}},$$

$$C_{\rm s}(r)=\big{<}{\bf S_{i}}\cdot{\bf S}_{{\bf i}+r{\bf e}_{x}}\big{>}\,,$$

$$C_{ijkl}=4\left(⟨({\mathbf{S}}_i\cdot {\mathbf{S}}_j)({\mathbf{S}}_k\cdot {\mathbf{S}}_l)⟩-⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_j{⟩}^2\right)$$

$$S_{\mathrm{VBC}}=\frac{1}{N_b}\sum _{k,l}\epsilon (k,l)C_{ijkl}$$

$$S_{\mathrm{VBC}}^{\ast }=\frac{1}{N_b}\sum _{k,l}^{\prime }\epsilon (k,l)C_{ijkl}$$

Full text

Critical colored-RVB states in the frustrated quantum Heisenberg model on the square lattice

Didier Poilblanc1*, Matthieu Mambrini1, Sylvain Capponi1

1 Laboratoire de Physique Théorique, C.N.R.S. and Université de Toulouse, 31062 Toulouse, France

• didier.poilblanc@gmail.com

Abstract

**We consider a family of SU(2)-symmetric Projected Entangled Paired States (PEPS) on the square lattice, defining colored-Resonating Valence Bond (RVB) states, to describe the quantum disordered phase of the $J_{1}-J_{2}$ frustrated Heisenberg model. For $J_{2}/J_{1}\sim 0.55$ we show the emergence of critical (algebraic) dimer-dimer correlations – typical of Rokhsar-Kivelson (RK) points of quantum dimer models on bipartite lattices – while, simultaneously, the spin-spin correlation length remains short. Our findings are consistent with a spin liquid or a weak Valence Bond Crystal in the neighborhood of an RK point. **

Contents

1. 1 Introduction: RVB and the frustrated Heisenberg model

2. 2 Numerical implementation

3. 2.1 iPEPS method

4. 2.2 Colored-RVB states

5. 2.3 CTMRG algorithm

6. 2.4 Optimization over the tensor parameters

7. 3 Results

8. 3.1 Energetics

9. 3.2 Dimer-dimer correlations

10. 3.3 Spin-spin correlations

11. 4 Conclusion and outlook

12. A CTMRG method

13. B VBC order parameters computed by ED

1 Introduction: RVB and the frustrated Heisenberg model

Resonant valence bond (RVB) states were first proposed by Anderson [1] to describe a possible spin liquid ground state of the $S=1/2$ antiferromagnetic Heisenberg model on the triangular lattice. Later on, it was also introduced as the parent Mott state of high-$T_{c}$ superconductors [2]. Soon after, along the same spirit, the Rokhsar Kivelson (RK) wavefunction [3] was defined as equal weight superposition of nearest neighbor (NN) dimer coverings, avoiding an explicit reference to the (hidden) spin degrees of freedom. It was shown that the RK wavefunction is a critical dimer liquid state [4] on the square lattice, in contrast to the case of non-bipartite kagome and triangular lattices [5, 6, 7, 8] on which a gapped ($\mathbb{Z}_{2}$) dimer liquid state is realized instead. Similarly, the NN RVB state, defined as an equal weight superposition of (non-orthogonal) NN singlet bond (also dubbed “dimer”) coverings, was shown to be also critical on the square lattice [9, 10] while several numerical work [11, 12, 13, 14] have demonstrated that their analogs on the kagome and triangular lattices are $\mathbb{Z}_{2}$ spin liquid states. Note that the (dimer) critical RK point is commonly unstable – ie towards dimerized phases [3, 15] or gapped dimer liquid phases [16] – upon slightly varying the model parameters. In fact, generically the RK point appears to be a multi-critical point with all sorts of nearby phases in which the critical correlations present at the RK point could be correct over a substantial intermediate range of energy and length scales. SU(2)-invariant spin models have also been engineered [17, 18] to mimic quantum dimer physics on the square lattice, with (critical) RVB ground state and Valence Bond Crystal (VBC) phases (spontaneously breaking translation symmetry), reflecting also the multi-critical nature of the RK point in SU(2)-symmetric systems.

Spin liquid behaviors are expected in two-dimensional (2D) frustrated quantum magnets where magnetic frustration prohibits magnetic ordering at zero temperature. Strong magnetic frustration is realized in the square lattice $J_{1}-J_{2}$ spin-1/2 Heisenberg model defined by summing over a 2D grid of lattice points $(i,j)$,

[TABLE]

and including both NN and next nearest neighbor (NNN) antiferromagnetic couplings $J_{1}$ (set to 1) and $J_{2}$, respectively. A paramagnetic quantum disordered (QD) region was suggested by early Lanczos Exact Diagonalizations (ED) extrapolations (including up to $N=36$ spins) in the range $J_{2}\in[0.34,0.68]$ [19], and similar results were announced later using ED up to $N=40$ [20]. However, until now, no agreement has been reached between several numerical approaches on the nature of the QD region – with proposals of VBC [21, 19, 22, 23, 24], (topological) gapped [25, 26] or gapless [27, 28, 29, 24] spin liquids. Interestingly, density matrix renormalization group (DMRG) approaches [30] with explicit implementation of SU(2) spin rotation symmetry [24] suggest that the QD region splits into a (critical) spin liquid phase (for $0.44<J_{2}<0.5$) and a plaquette VBC phase (for $0.5<J_{2}<0.61$). Recently, DMRG simulations of Wang and Sandvik using level spectroscopy [31] also indicate that the QD is formed by a gapless spin liquid phase (for $0.46<J_{2}<0.52$) and a VBC (for $0.52<J_{2}<0.62$). In contrast, other recent computations using U(1)-symmetric (infinite size) Projected Entangled Pair States (PEPS) [32] suggest a columnar VBC (for $0.53<J_{2}<0.61$) separated from the conventional Néel phase by a deconfined critical point [33], in qualitative agreement with a previous finite size PEPS computation [34].

Despite such recent progress, the exact nature of the QD phase remains still unclear. In this paper we aim to investigate further the QD phase in the region around $J_{2}=0.55$ introducing simple PEPS Ansätze which are specially designed to describe SU(2)-invariant states with full space group symmetry. In Sec. 2 we quickly review the iPEPS method used, further details being provided in Appendix A. Results on variational energy and correlation functions are analyzed in Sec. 3 and complementary ED results are provided in Appendix B, strongly suggesting that RK physics with long-range dimer correlations emerges. Finally, further discussions and conclusions are given in the last section 4.

2 Numerical implementation

2.1 iPEPS method

Tensor networks [35, 36, 37, 38] have recently emerged as a state-of-the-art numerical tool to tackle correlated lattice models. Among them, 2D PEPS [39, 40] are variational Ansätze constructed from local site tensors carrying the physical degrees of freedom (of dimension $2$ for spin-$\frac{1}{2}$ systems) and $z$ “virtual” bonds ($z$ is the lattice coordination number, $z=4$ for the square lattice) of arbitrary dimension $D$ (see Appendix (A)). Interestingly, local (gauge) or global (physical) symmetries can be implemented in PEPS [41, 42, 43, 44, 45, 46, 47, 48]. In the infinite-PEPS (iPEPS) method [49], one works directly in the thermodynamic limit by approximating the (infinite) space around a small $M$-site cluster by an effective “environment” (here $M=2\times 2=4$). One of the most accurate computation of the environment is based on a Tensor Renormalization Group (TRG) scheme involving Corner Transfer Matrices (CTM) [50, 51, 52, 53]. Unrestricted energy minimization can be performed [54, 55] using a simple update [56, 57] or a full update [58] of the environment. Recently, a new variational optimization scheme using a Conjugate Gradient (CG) algorithm has been tested on the non-frustrated [59, 60] and on the above spin-1/2 $J_{1}-J_{2}$ Heisenberg model [61, 62].

2.2 Colored-RVB states

We wish here to refine and extend the previous iPEPS study of Ref. [62] dealing with the spin-1/2 frustrated $J_{1}-J_{2}$ Heisenberg model on the square lattice. While Ref. [62] focused mainly on $J_{2}=0.5$ – pointing towards a gapless spin liquid – we focus here on slightly larger $J_{2}\sim 0.55$ where, we shall argue, a new behavior occurs. For this purpose, we shall consider the same families of translationally invariant fully symmetric PEPS involving a linear combination of a finite number $\cal D$ of single site tensors,

[TABLE]

The tensors $t_{\alpha}$, obtained from a complete classification of symmetric site tensors on the square lattice [48], are fully invariant under SU(2) spin rotations and under all operations of the $C_{4v}$ point group (90-degree rotations and reflections). These local symmetry properties of the site tensors guarantee that the PEPS itself is a fully symmetric wavefunction under all the global symmetry operations leaving the Hamiltonian invariant. The bond virtual space of dimension $D=2{\cal N}+1$ is of the form

[TABLE]

corresponding to $\cal N$ possible “colors” of spin-$\frac{1}{2}$ and a spin-[math] (singlet). In the following we shall consider the three PEPS families associated to one, two and three colors of the spin-1/2 degree of freedom i.e. namely to $V=\frac{1}{2}\oplus 0$ (${\cal N}=1,D=3$), $V=\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ (${\cal N}=2,D=5$), and $V=\frac{1}{2}\oplus\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ (${\cal N}=3,D=7$). Each of these PEPS family is spanned by a small number $\cal D$ of linearly independent tensors, ${\cal D}=2$, 10 and 30 respectively, given in Refs. [48, 62].

It is important to notice that these three PEPS families are not separated from each other but rather embedded into one another. The smallest one (spanned by 2 independent $D=3$ site tensors) can be viewed as a manifold of generalized RVB states which include, when expended in terms of valence bond (VB) configurations, singlet bonds extending beyond NN sites (in contrast to the original NN RVB state [1]). Its corresponding phase diagram contains a RK dimer liquid phase and a (topological) spin liquid phase [63]. The ${\cal N}\geq 2$ PEPS family can be viewed as $\cal N$-color (generalized) RVB states where singlet valence bonds (VB) carry now a color index, ranging from 1 to $\cal N$, and where the VB amplitudes depend on the coloring pattern. In that new language, it becomes obvious that the ${\cal N}=2$ PEPS family includes two (disjoint) copies of the manifold of single-color RVB states. Similarly, our largest ${\cal N}=3$ PEPS family – defining the manifold of 3-color RVB states – contains three (disjoint) sub-manifolds of single-color RVB states (of different colors) and three (disjoint) sub-manifolds of 2-color RVB states (with different pairs of colors). These features are summarized in Fig. 1.

All site tensors contain either one or three spin-0 virtual legs (leading to a $\mathbb{Z}_{2}$ gauge symmetry [63] associated to the odd parity of this number of legs). Note that, if one restricts to only the subset of tensors with a single spin-1/2 leg and three spin-0 legs, the corresponding PEPS is the usual NN RVB state (the VB amplitudes do not depend on the coloring pattern in that case). Longer range singlets are created by “teleportation” [28, 63] introduced by any of the site tensors containing three spin-1/2 (of any color) and one spin-0 on the virtual legs. Hence, the most general 3-color ($D=7$) RVB Ansatz can be viewed as a resonant state of colored VB coverings of the type drawn in Fig. 2. The VB amplitude depends both on the VB covering and on the coloring pattern in a complex way set up by the tensor coefficients $c_{\alpha}$ entering Eq. 2.

2.3 CTMRG algorithm

For a given PEPS realization (i.e. defined by a particular set $\{c_{\alpha}\}$ of coefficients in Eq. 2) the corresponding energy $E\left[\{c_{\alpha}\}\right]$ (in the thermodynamic limit) is computed by a CTMRG method which takes advantage of the point group symmetry of the lattice (see Appendix (A) for details). Note that, although the site tensor is fully SU(2)-invariant, the CTMRG procedure of Ref. [62] (used to contract the infinite tensor network outside a $2\times 2$ active region) was generically converging to a fixed-point environment exhibiting a small finite staggered magnetization i.e. spontaneously breaking SU(2) spin-rotation symmetry, at least for $J_{2}=0.5$ which was extensively studied. Although, this effect is spurious (the data are consistent with a vanishing staggered magnetization in the limit of infinite environment dimension, $\chi\rightarrow\infty$), it complicates the analysis of the data. Hence, we have improved the CTMRG procedure in order to keep the full SU(2) symmetry at all stage and for all $\chi$ (despite numerical rounding errors) so that the fixed-point solution of the CTMRG indeed corresponds to a fully symmetric QD state.

2.4 Optimization over the tensor parameters

In summary, the CTMRG algorithm enables to compute the energy density $E\left[\{c_{\alpha}\}\right]$, in the thermodynamic limit, for a given choice of (i) the bond dimension $D$, (ii) the associated $\cal D$ tensor coefficients $c_{\alpha}$ and (iii) the environment dimension $\chi$. We use a brute force (Conjugate Gradient) optimization upon the set of coefficients $c_{\alpha}$ to obtain the absolute minimum of the variational energy at given $D$ and $\chi$. This requires to numerically compute each component of the local gradient vector $\vec{\cal G}$ of the energy $E\left[\{c_{\alpha}\}\right]$ by finite differentiation,

[TABLE]

where, in the set of parameter $\{c_{\alpha}\}_{{}_{\beta}}$, only $c_{\beta}$ has been incremented to $c_{\beta}+\delta$, $\delta/c_{\beta}$ corresponding typically to a relative change of less than $1\%$. Note that in the calculation of $E\left[\{c_{\alpha}\}_{{}_{\beta}}\right]$ it is crucial to take into account the change of the environment by computing the new CTMRG fixed point. Interestingly, we observe that the color-exchange symmetry is broken in the optimal PEPS. The optimization is performed for each choice of ${\cal N}=1,2,3$ and up to some maximum value of $\chi$, $\chi_{\rm opt}(D)$.

Thanks to the refinements of the iPEPS technique mentioned above, we have obtained accurate results for $J_{2}\sim 0.55$ detailed below.

3 Results

3.1 Energetics

The variational energies (per site) at $J_{2}=0.55$ are shown in Fig. 3(a) as a function of the inverse of the environment dimension $\chi$. Note that the local tensors are fully optimized up to a maximum bond dimension $\chi_{\rm opt}=12D^{2}=108$ for $D=3$, $\chi_{\rm opt}=4D^{2}=100$ for $D=5$ and $\chi_{\rm opt}=2D^{2}=98$ for $D=7$. Then, energies are computed for larger environment dimensions $\chi>\chi_{\rm opt}$ using the fixed optimized tensors obtained at $\chi=\chi_{\rm opt}$. Linear fits can be performed in $1/\chi$ to provide the $\chi\rightarrow\infty$ true variational energies. Note that our energy $-0.4842$ for $D=7$ is quite close to the value $-0.4856(1)$ obtained using finite $D=9$ PEPS cluster update  [34] and finite $D=8$ PEPS Variational Monte Carlo (VMC) [64, 65]. Moreover, the $D\rightarrow\infty$ extrapolation shown in Fig. 3(b), using either a Taylor series or a power-law in $1/D$, gives $-0.4894(5)$, significantly below the DMRG [24] and the VMC [29] estimates (reported in Fig. 3(b) for convenience). Alternatively, a power-law extrapolation (almost linear) w.r.t. $1/{\cal N}$ gives a slightly lower energy $-0.4909$. This gives us some confidence that the series of colored-RVB states, as defined by the PEPS construction, provides a faithful representation of the low-energy physics of the $J_{1}-J_{2}$ model at frustration $J_{2}\sim 0.55$.

3.2 Dimer-dimer correlations

After optimizing our symmetric PEPS Ansatz w.r.t. the coefficients of the site tensor, correlation functions can be computed using arbitrarily long (let’s say horizontal) one dimensional strips bounded by environment tensors [62] (which depend on $\chi$) as depicted in Fig.7 (e). Let us first define the connected dimer-dimer correlations,

[TABLE]

where ${\bf x}=(i,j)$ is some arbitrary lattice site, dimer operators $D_{\bf x}={\bf S_{x}}\cdot{\bf S}_{{\bf x}+{\bf e}_{x}}$ are oriented along the horizontal ${\bf e}_{x}=(1,0)$ direction, and the expectation values are taken in the optimized PEPS. Note that, the PEPS being invariant by lattice translation, the dimer density \big{<}D_{\bf x}\big{>} does not depend in fact on the position ${\bf x}$. Also, although we are using a strip geometry, the local tensor (and the corresponding environment tensor $T$) has been optimized for the fully rotationally invariant (infinite) lattice.

The dimer correlations are plotted in Fig. 4(a) for $D=7$ and $J_{2}=0.55$ and several values of $\chi$, in semi-log scale to reveal the long-distance exponential decay. From a linear fit, one can extract the corresponding dimer correlation length $\xi_{d}(\chi)$. The latter is plotted in Fig. 5(a) as a function of $\chi$ and, in Fig. 5(b), versus $\chi/D^{2}$ which seems to be the natural rescaled variable to compare the behaviors of the 3 different families $D=3,5,7$. For all cases, we observe a clear linear dependence with $\chi$,

[TABLE]

suggesting a divergence of the correlation length, or at least saturation to a very large value beyond reach. It is interesting to notice that, once plotted in terms of the rescaled variable $\chi/D^{2}$, the slope $a_{D}$ of the linear increase is quite similar for $D=5$ and $D=7$, suggesting a robust feature of the correlations. In Fig. 5(c) we compare the $D=7$ correlation length at different $J_{2}$ values, showing a more pronounced increase at $J_{2}=0.55$.

Whenever the correlation length $\xi_{d}(\chi)$ diverges (or becomes very large), one expect to see power-law behaviors in the correlation functions,

[TABLE]

in the range of distance $1<r<\xi_{d}$. The exponent can be written as $\alpha_{\rm d}=1+\eta_{\rm d}$ where $\eta_{\rm d}$ defined e.g. in Ref. [66] is the anomalous dimension. Since the correlation length remains moderate for $D=7$, we have (i) first extrapolated the data in the $\chi\rightarrow\infty$ limit in Fig. 4(b) (using a power-law fit in $1/\chi$) for a few distances $r$ and (ii) fitted these extrapolated values to obtain the exponent $\alpha_{d}\sim 1.08$ via a power-law fit in Fig. 4(a). The smallness of the anomalous exponent $\eta_{d}\simeq 0.08$ reveals very long-range dimer correlations at $J_{2}=0.55$, in contrast to $J_{2}=0.5$ studied in Ref. [62]. Interestingly, quite similar behaviors are found for $D=3,5$ and $7$ as shown in Fig. 4(c) where the dimer correlations are plotted, using a log-log scale, as a function of the rescaled distance ${\tilde{r}}=r/\xi_{d}$, for the largest attainable environment dimension $\chi$. Linear fits for $\tilde{r}<1$ provide similar values for the exponent $\alpha_{d}$, between $1.15$ and $1.25$, in agreement (within error bars) with the previous analysis. It is also interesting to notice that these values are quite close to the value $\alpha_{\rm d}\simeq 1.16$ reported for the NN RVB state [9] and agree with recent DMRG simulations [24].

3.3 Spin-spin correlations

Finally, we have computed the spin-spin correlations (e.g. along the ${\bf e}_{x}$ horizontal direction),

[TABLE]

using the same strip geometry of an (infinite) chain of sites bounded by environment tensors $\cal T$ on the edges (Fig.7 (d)). In the original RVB picture [1] spins are correlated only at short distance via NN singlet pairing. Our results obtained for $J_{2}=0.55$, $D=7$ and several $\chi$ values up to $\chi=9D^{2}=441$ are shown in Fig. 6. Linear fits of the long-distance correlations (plotted in log-log scales) enable to estimate accurately the spin correlation length. The inset shows that the latter remains quite small, typically less than 2 lattice spacings, even for the largest $\chi$ at hand. The same is also true for $D=3$ and $D=5$. However, as shown in Fig. 4(a), the spin correlations are much stronger than the dimer correlations at short distance. This is consistent with the RVB picture where strong (resonating) singlet bonds are formed between NN sites.

4 Conclusion and outlook

The main findings of this iPEPS study at $J_{2}=0.55$ and $J_{2}=0.575$ are the following: i) a simple SU(2)-symmetric PEPS based on a single site tensor provides a very good variational energy; ii) its dimer correlations exhibit slow algebraic decay up to long distance; iii) its spin correlations are short range. Properties ii) and iii) are characteristic of RK physics found e.g. in the NN RVB spin liquid on any bipartite lattice.

The critical RK point is known to be unstable to small Hamiltonian perturbations [7, 15, 16] breaking the lattice bipartiteness. However, investigation of classical dimer models at finite temperature [8] indicates that criticality and nonbipartiteness are compatible. In fact, our $D=3$ PEPS, which realizes exactly an extended RVB state with (inter-sublattice) longer-range singlets, is known to possess an extended RK dimer liquid phase [63]. Also, it is likely that regions of (truly critical) RK phases exist also within our $D=5$ and $D=7$ PEPS manifolds. Although one cannot prove that such a RK dimer liquid is realized in the $J_{1}-J_{2}$ spin-1/2 Heisenberg antiferromagnet, it is known that the critical RVB state is the ground state of a family of SU(2)-symmetric local spin-1/2 models with frustrating interactions [17, 18]. In any case, the critical dimer correlations could survive in nearby phases of some RK point over a substantial intermediate range of distances. In that case, (at least) two scenario (probably beyond our current PEPS description) may apply; First, it may well be that the dimer correlation length saturates to a (very) large value leading to a (gapped) spin liquid with $\xi_{d}\gg\xi_{s}$. A second possibility is that the system would spontaneously break translation symmetry and develops a type of (very weak) VBC ordering (dimerization, plaquette formation, …) as suggested by large-N theories [67], series expansions [68, 69] or numerical work [21, 22, 24, 32, 31]. In fact, Lanczos ED of small clusters (see Appendix (B) for details) suggests that the tendency to realize a VBC is maximum at $J_{2}\simeq 0.55$, although the VBC order parameter should be quite small, and probably very hard to detect directly. Note that a gapped spin liquid with $\xi_{d}\gg\xi_{s}$ could be alternatively seen as a melted VBC.

Acknowledgments

We acknowledge useful discussions with Fabien Alet and thank Juraj Hasik and Laurens Vanderstraeten for private communications. D.P. also acknowledges inspiring conversations with Federico Becca, Paul Fendley, Zheng-Cheng Gu, Steve Kivelson, and Wen-Yuan Liu. S. C. thanks Alex Wietek for making some of the ED simulations possible using MPI technique [70].

Funding information

This project is supported by the TNSTRONG ANR-16-CE30-0025 and TNTOP ANR-18-CE30-0026-01 grants awarded by the French Research Council. This work was granted access to the HPC resources of CALMIP supercomputing center (under the allocations P1231 and P0677) and GENCI (Grant number A0050500225).

Appendix A CTMRG method

In this appendix, we provide a brief and self-contained description of the CMTRG method used to characterize the properties of the iPEPS states considered in this paper. We focus on the renormalization procedure aiming at deriving converged environment tensors (corner and edge tensors) at the thermodynamic limit that can be further used to compute states properties such as energy or correlation functions. The discussion is restricted to the case of a fully symmetric tensor (i.e. transforming according to the $A_{1}$ representation of $C_{4v}$) in the context of a translationally invariant Ansatz (the same tensor is used on every site of the square lattice). For a more general presentation, one can refer to the appendix A of reference [53].

Bilayer tensors. In the infinite-PEPS (iPEPS) method [49], one considers a PEPS Ansatz |\Psi\big{>} directly in the thermodynamic limit. The PEPS is an infinite two-dimensional tensor network defined by a single site tensor, and its normalization \big{<}\Psi|\Psi\big{>} is then a bilayer tensor network which can be re-expressed as a tensor network of site rank-4 bilayer tensors (of bond dimension $D^{2}$). The bilayer tensor is represented in Fig. 7(a) and possesses full invariance under spin rotation and point group symmetry operations.

Observables. Computation of observables (like energy or correlations) also requires the bilayer tensor network which is approximately contracted over the (infinite) space surrounded a small $M$-site cluster. This approximate contraction then leads to an effective “environment” of this small region. For the energy one needs a $M=2\times 2=4$ site cluster (fitting the interaction on both NN and diagonal bonds, see Fig. 7(c)) and, for the correlations at distance $r$, a one-dimensional $r$-site segment connecting two operators at its two ends (see Fig. 7(d) and (e)).

Renormalization procedure. The computation of the environment is based on a Corner Transfer Matrix Renormalization Group (CTMRG) [50, 51, 52, 53] scheme schematically represented in Fig. 7(b).

The environment involves a $\chi\times\chi$ corner transfer matrix $\cal C$ and a rank-3 boundary $\chi\times\chi\times D^{2}$. In practice $\chi=kD^{2}$ with $k$ integer. Before describing the several steps of the CTRMG algorithm let us remark that, thanks to the $A_{1}$ symmetry of the site tensor, several important simplifications occurs in the procedure. First of all, the four corner matrices as well as the four edge tensors are degenerate, so that a single $({\cal C},{\cal T})$ couple is needed. Furthermore, $\cal C$ is a real symmetric matrix. Hence it can be reduced by diagonalization (instead of a singular value decomposition) and only one isometry $\cal U$ has to be considered.

Initialization step. Corner matrix $\cal C$ and edge tensor $\cal T$ are initialized in a similar way as bilayer site tensor is constructed from the site tensor (Fig. 7(a)). In addition to the physical index, one (resp. two) virtual bonds are contracted between the two layers. 2. 2.

Corner renormalization. The new corner $\tilde{\cal C}$ is obtained in two steps. Starting with a $\cal TCT$ corner, one adds a bilayer site tensor $\cal A$ (see yellow square on Fig. 7(b)). The resulting (real) symmetric $\chi D^{2}\times\chi D^{2}$ matrix is diagonalized and an $\chi D^{2}\times\chi$ isometry $\cal U$ is constructed by keeping only (at most) the $\chi$ largest weights. Special care is taken to preserve the SU(2) spin-rotation symmetry in the truncation by keeping the SU(2) multiplet structure appearing in the corner matrix spectrum. 3. 3.

Edge Renormalization By adding a bilayer site tensor $\cal A$ to the edge tensor $\cal T$ and contracting with the isometry $\cal U$, the renormalized $\chi\times\chi\times D^{2}$ edge tensor $\tilde{\cal T}$ is obtained (see yellow rectangle on Fig. 7(b)).

Steps 2. and 3. are then repeated until a fixed point for ${\cal C}$ is reached. Note that the complexity of step 2. is $\chi^{3}D^{2}+\chi^{3}D^{4}+\chi^{2}D^{8}=k^{3}D^{8}+k^{3}D^{10}+k^{2}D^{12}\sim k^{2}D^{12}$ for the untruncated corner matrix computation and $\left(\chi D^{2}\right)^{3}=k^{3}D^{12}$ for the diagonalization. The cost of step 3. is $\chi^{3}D^{4}+\chi^{3}D^{6}+\chi^{2}D^{8}=k^{3}D^{10}+k^{2}D^{12}+k^{3}D^{12}\sim k^{3}D^{12}$. As a result, the algorithmic complexity is $D^{12}$ for large $D$.

Appendix B VBC order parameters computed by ED

In order to investigate if the ground-state is a VBC that breaks lattice symmetries, we have computed the dimer-dimer correlation function:

[TABLE]

on various finite-size tori of $N$ sites. Following Ref. [22], we can then compute various structure factors, and in particular

[TABLE]

where the summation is over $N_{b}$ parallel bonds $(kl)$ with respect to the reference bond $(ij)$ and $\varepsilon(k,l)=\pm 1$ depending on the sublattice, see Fig. 8(c). It can be shown that $S_{\mathrm{VBC}}$ is finite both for a VBC with columnar or plaquette order [22].

In Fig. 8(a), we plot the behavior of $S_{\mathrm{VBC}}$ vs $J_{2}$ for different finite-size clusters. In order to remove short distance data from this order parameter, we have also considered a slightly different definition

[TABLE]

where the summation does not include the nearest six bonds, see Fig. 8(b,c). Note that, compared to the ground-state energy calculations, we only computed VBC order parameter on clusters that are compatible with plaquette or columnar order, i.e. contain $(\pi,0)$ and $(0,\pi)$ in their Brillouin zone. Quite interestingly, both VBC order parameters are maximal around $J_{2}\simeq 0.55$, which is the optimal value found in DMRG [24], and then have a sudden drop beyond $J_{2}=0.6$, which is presumably of first-order character.

Finite-size scaling analysis is shown in Fig. 9 for VBC order parameters at various $J_{2}$ values. Reliable extrapolation is not possible, but given the data points and their curvature vs $1/N$, data are compatible with a vanishing VBC order parameter for $J_{2}=0.4$, $0.5$ or $0.6$, but weak long-range VBC order could be stabilized around $J_{2}=0.55$.