Quantum critical phase with infinite projected entangled paired states

TL;DR

This paper introduces a class of SU(2)-invariant PEPS that effectively models a quantum critical phase in the frustrated J1-J2 Heisenberg model, showing no magnetic order and diverging correlation lengths.

Contribution

It identifies a specific SU(2)-invariant PEPS family that captures quantum critical behavior in the J1-J2 model, advancing understanding of quantum phase transitions.

Findings

No magnetic order observed in the PEPS for N≥3

Diverging correlation lengths indicating criticality

Effective variational energies at J2=0.5

Abstract

A classification of SU(2)-invariant Projected Entangled Paired States (PEPS) on the square lattice, based on a unique site tensor, has been recently introduced by Mambrini et al.~\cite{Mambrini2016}. It is not clear whether such SU(2)-invariant PEPS can either i) exhibit long-range magnetic order (like in the N\'eel phase) or ii) describe a genuine quantum critical point (QCP) or quantum critical phase (QCPh) separating two ordered phases. Here, we identify a specific family of SU(2)-invariant PEPS of the classification which provides excellent variational energies for the $J_1-J_2$ frustrated Heisenberg model, especially at $J_2=0.5$, corresponding to the approximate location of the QCP or QCPh separating the N\'eel phase from a dimerized phase. The PEPS are build from virtual states belonging to the $\frac{1}{2}^{\otimes N} \oplus 0$ SU(2)-representation, i.e. with $N$ "colors" of virtual \hbox{spin-$\frac{1}{2}$}. Using a full update infinite-PEPS approach directly in the thermodynamic limit, based on the Corner Transfer Matrix renormalization algorithm supplemented by a Conjugate Gradient optimization scheme, we provide evidence of i) the absence of magnetic order and of ii) diverging correlation lengths (i.e. showing no sign of saturation with increasing environment dimension) in both the singlet and triplet channels, when the number of colors $N\ge 3$. We argue that such a PEPS gives a qualitative description of the QCP or QCPh of the $J_1-J_2$ model.

Tables

Table 1. Table 1: [Color online] List of all virtual spaces V 𝑉 V of bond dimension D ≤ 7 𝐷 7 D\leq 7 for which V ⊗ 4 superscript 𝑉 tensor-product absent 4 V^{\otimes 4} can be projected onto a physical spin 1/2. The ones considered here are indicated by (green) marks. Classes with higher spins give poorer variational energies than the lower spin ones of same total bond dimension D 𝐷 D .

$D$	3		4		5		6	7
$V$	✓	$\frac{1}{2}\oplus 0$	✓	$\frac{1}{2}\oplus 0\oplus 0$		$\frac{1}{2}\oplus 0\oplus 0\oplus 0$	$1\oplus \frac{1}{2}\oplus 0$	✓	$\frac{1}{2}\oplus \frac{1}{2}\oplus \frac{1}{2}\oplus 0$
					✓	$\frac{1}{2}\oplus \frac{1}{2}\oplus 0$			$1\oplus \frac{1}{2}\oplus 0\oplus 0$
					✓	$1\oplus \frac{1}{2}$			$\frac{3}{2}\oplus \frac{1}{2}\oplus 0$
						$\frac{3}{2}\oplus 0$			$\frac{3}{2}\oplus 1$
									$2\oplus \frac{1}{2}$
									$\frac{5}{2}\oplus 0$

Table 2. Table 2: Comparison between our D = 7 𝐷 7 D=7 iPEPS results ( χ → ∞ → 𝜒 \chi\rightarrow\infty extrapolations) and the best estimates in the literature, for J 2 = 0 subscript 𝐽 2 0 J_{2}=0 , J 2 = 0.5 subscript 𝐽 2 0.5 J_{2}=0.5 and J 2 = 0.55 subscript 𝐽 2 0.55 J_{2}=0.55  : J 2 = 0 subscript 𝐽 2 0 J_{2}=0 results are obtained by QMC  Sandvik and Evertz ( 2010 ); Sandvik ( 2010 ) . At finite J 2 subscript 𝐽 2 J_{2} , we quote energies obtained by extrapolations to the thermodynamic limit using DMRG  Gong  et al. ( 2014 ) , VMC  Hu  et al. ( 2013 ) and finite-size D = 9 𝐷 9 D=9 PEPS  Wang  et al. ( 2016 ) . Note that the D = 7 𝐷 7 D=7 iPEPS energies are only upper bounds of the true variational energies (see text).

J	0	0.5	0.55
QMC	$-0.66944$
DMRG		$-0.4968$	$-0.4863$
VMC		$-0.4970\left(5\right)$	$-0.4870\left(5\right)$
$D=9$ PEPS		$-0.4958\left(3\right)$	$-0.4857\left(2\right)$
$D=7$ iPEPS	$-0.6677$	$-0.4950$	-0.4830

Equations

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

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

$$C_{\rm d}^{(\rm T)}(d)=\big{<}D_{\bf i}^{y}D_{{\bf i}+d{\bf e}_{x}}^{y}\big{>}_{0}-\big{<}D_{\bf i}^{y}\big{>}_{0}\big{<}D_{{\bf i}+d{\bf e}_{x}}^{y}\big{>}_{0}\,,$$

$$C_{\mathrm{s}}(d)$$

$$C_{\mathrm{d}}(d)$$

$${\xi }_D(\chi )\simeq f_D\chi /D^2,$$

$$C_{\lambda }(d)\sim C_0\exp (-d/{\xi }_D(\chi )),$$

$${\mathbf{S}}_{\mathbf{i}}=S_{\mathbf{i}}^{\parallel }\mathbf{n}+{\mathbf{S}}_{\mathbf{i}}^{\perp },$$

$$C_{\mathrm{s}}^{\parallel }(d)$$

$$C_{\mathrm{s}}^{\perp }(d)$$

Full text

I Introduction

Low-dimensional quantum magnets offer a rich zoo of phases breaking a discrete (like point group or lattice) or a continuous (like spin rotation) symmetry. Often, such phases are separated by Quantum Critical Points (QCP), as described within the usual Ginsburg-Landau (GL) framework. Interestingly, it has been proposed that some QCP may not be described by the GL paradigm Senthil et al. (2004); Shao et al. (2016). A celebrated quantum spin model is the frustrated spin-$\frac{1}{2}$ Heisenberg model on the two-dimensional (2D) square lattice involving competition between nearest neighbor (NN) and next-nearest neighbor (NNN) antiferromagnetic (AF) couplings, $J_{1}$ and $J_{2}$ respectively. Setting $J_{1}=1$, $J_{2}$ controls the amount of frustration which is maximum (classically) at $J_{2}=0.5$. Large-scale Quantum Monte Carlo (QMC) simulations Reger and Young (1988); Sandvik and Evertz (2010); Sandvik (2010) has shown that the ground state (GS) of the unfrustrated ($J_{2}=0$) Heisenberg model exhibits long range (LR) AF order. In the thermodynamic limit, the (global) spin-rotational SU(2) symmetry is spontaneously broken and the GS acquires a finite local staggered magnetization. When $J_{2}$ is turned on, the order parameter is gradually suppressed and a quantum phase transition to a Quantum Disordered (QD) phase Chandra and Douçot (1988); Dagotto and Moreo (1989a); Sachdev and Read (1991); Mila et al. (1991) – such as a dimer Read and Sachdev (1989); Dagotto and Moreo (1989b); Poilblanc et al. (1991); H.J. Schulz et al. (1996) or a plaquette Zhitomirsky and Ueda (1996); Mambrini et al. (2006) Valence Bond Crystal (VBC) – takes place (see Fig. 1). It was also argued that magnetic frustration could stabilize spin liquids (with no symmetry breaking), such as the Resonating Valence Bond (RVB) states Anderson (1973) showing algebraic (short range) VBC correlations on the square (Kagome) lattice Albuquerque and Alet (2010); Tang et al. (2011); Poilblanc et al. (2012); Schuch et al. (2012).

Recently, tremendous progress have been made in tensor network techniques Cirac and Verstraete (2009); Cirac (2010); Orús (2014a); Schuch (2013); Orús (2014b), aiming to go beyond Density Matrix Renormalization Group (DMRG) methods White (1992) in 2D. More specifically, Projected Entangled Pair States (PEPS) Pérez-García et al. (2008) are variational ansätze constructed from a few local tensors, located on $M$ non-equivalent sites, and characterized by (i) one bond carrying the physical degrees of freedom (of dimension $2$ for spin-$\frac{1}{2}$ systems) and (ii) $z$ “virtual” bonds ($z$ is the lattice coordination number, $z=4$ for the square lattice) of arbitrary dimension $D$ as shown in Fig. 2(a). Interestingly, any local (gauge) or global (physical) symmetry can be implemented in PEPS Pérez-García et al. (2010); Schuch et al. (2010); Singh et al. (2010); Singh and Vidal (2012, 2013); Weichselbaum (2012); Jiang and Ran (2015); Haegeman et al. (2015); Mambrini et al. (2016). Also, a simple bulk-edge (holographic) correspondence provides a remarkable tool to investigate the properties of edge states Cirac et al. (2011); Yang et al. (2014). Many remarkable states of matter such as trivial paramagnets Jian and Zaletel (2016), topological Schuch et al. (2012); Poilblanc et al. (2012); Poilblanc and Schuch (2013) or algebraic Poilblanc et al. (2012) RVB spin liquids, loop spin liquids Li et al. (2014), superfluids Poilblanc et al. (2013) or unconventional correlated superconductors Poilblanc et al. (2014) have simple representations in terms of PEPS. Numerical calculations with PEPS do not require to compute the wave function coefficients (which, conceptually, are given by contracting the tensor network over all virtual links) but, rather, make use of transfer matrices Haegeman and Verstraete based on “double-layer” tensors (see Fig. 2(b-e)). In the infinite-PEPS (iPEPS) method Jordan et al. (2008), one works directly in the thermodynamic limit by approximating the (infinite) space around a small $M$-site cluster by an effective “environment” (see Fig. 2(f)). One of the most accurate computation of the environment is based on a Renormalization Group scheme involving Corner Transfer Matrices (CTMRG) Nishino and Okunishi (1996); Nishino et al. (2001); Orús and Vidal (2009); Orús (2012) as shown in Fig. 2(g,h). Unrestricted energy minimization over the $MdD^{z}$ tensor coefficients can be performed using Time Evolution Block Decimation (TEBD) Vidal (2007); Orús and Vidal (2008) which has to be combined with a simple update Vidal (2003); Jiang et al. (2008) or a full update Phien et al. (2015) of the environment. A (finite) PEPS method using a $2\times 2$ cluster update supplemented by a finite size extrapolation has also been introduced Wang et al. (2016). Recently, a new optimization scheme using a Conjugate Gradient (CG) algorithm has been tested on the non-frustrated Corboz (2016); Vanderstraeten et al. (2016) and frustrated Liu et al. Heisenberg model, with iPEPS or finite PEPS, respectively.

The entanglement entropy (i.e. the quantity measuring the amount of entanglement in a bi-partitioned system) in a spontaneously-broken state exhibits anomalous additive logarithmic corrections Kallin et al. (2011); Metlitski and Grover (2011); Ju et al. (2012) to the area law (i.e. the linear scaling of the entropy with the length of the cut). When the staggered magnetization $m_{\rm stag}\rightarrow 0$, at the QCP, the violation of the area law is expected to be even more severe. This means that a good description of the QCP, or even of the Néel state, in terms of a PEPS (which strictly fulfills the area law for any finite $D$) is particularly challenging. A very simple ($D=3$) PEPS ansatz for the Néel state on the square lattice was first proposed in terms of a (one-parameter) spinon-doped RVB phase Poilblanc (2014). Also, finite size PEPS Lubasch et al. (2014) or, more recently, state-of-the-art iPEPS calculations involving a Conjugate Gradient (CG) minimization algorithm Corboz (2016); Vanderstraeten et al. (2016) came up with very accurate energy for the Néel GS of the 2D Heisenberg model. However, the phase diagram of the $J_{1}-J_{2}$ model is still heavily debated. No agreement has been reached between several numerical approaches, neither on the nature of the QD region – with proposals of VBC H.J. Schulz et al. (1996); Mambrini et al. (2006); Capriotti and Sorella (2000); Gong et al. (2014), (topological) gapped Jiang et al. (2012) or gapless Capriotti et al. (2001); Wang et al. (2013); Hu et al. (2013); Gong et al. (2014); Morita et al. (2015) spin liquids – nor on the location $J_{2}=J_{2c}$ of the phase transition. While early Exact Diagonalisations (ED) extrapolations H.J. Schulz et al. (1996) were bracketing $J_{2c}\in[0.34,0.6]$, DMRG studies Jiang et al. (2012); Gong et al. (2014) suggested $J_{2c}\simeq 0.41-0.44$, while Variational Monte Carlo (VMC) studies Hu et al. (2013); Morita et al. (2015) give $J_{\rm 2c}\simeq 0.48-0.5$ and finite-size (cluster update) PEPS computations Wang et al. (2016) $J_{2c}\simeq 0.572(5)$. Recently, Wang and Sandvik Wang and Sandvik (2017) argued for a quantum critical phase (QCPh) centered around 0.5. In all these approaches (except ED), the spin rotational SU(2) symmetry is explicitly broken in the Néel phase. However, there is no obstruction principle to construct accurate SU(2)-symmetric wave functions exhibiting long range AF order Kaneko et al. (2016). Since such states may be characterized by a large entanglement, it is unclear whether it can be realized with low-$D$ symmetric PEPS. Also, whether SU(2)-symmetric PEPS have the potential to describe zero-temperature QCP or QCPh – in the same way as one-dimensional (1D) Matrix Product States (MPS) can describe critical 1D systems Tagliacozzo et al. (2008); Pollmann et al. (2009); Pirvu et al. (2012) – is still unclear Liu et al. (2010). Though, it is known that non-trivial criticality can be captured by PEPS, even at finite $D$ Nishino et al. (2001); Verstraete et al. (2006).

Motivated by the above conceptual and practical issues, we have re-visited the $J_{1}-J_{2}$ model using some new PEPS developments, based on a general scheme to construct SU(2)-symmetric PEPS using computer-assisted algebra Mambrini et al. (2016). This enables us to introduce key features in the full-update iPEPS scheme : (i) Full translational and rotational invariance is enforced by using a unique SU(2)-invariant tensor on every lattice site; (ii) Full optimization of the (few) tensor coefficients is accomplished via a CG method; (iii) Careful scaling with environment dimension $\chi$ is performed in order to address the $\chi\rightarrow\infty$ limit. Using this procedure, we have identified a specific (low-dimensional) family of SU(2)-symmetric PEPS which provides excellent variational energies for the $J_{1}-J_{2}$ frustrated Heisenberg model, especially at $J_{2}=0.5$, i.e. close to the (unknown) QCP or QCPh of this model. We show evidence that these (optimized) PEPS do not exhibit long range AF order. We also find that, above bond dimension $D=7$, the PEPS (optimized for $J_{2}=0.5$) exhibits diverging spin-spin and dimer-dimer correlation lengths, i.e. showing no sign of saturation up to large environment dimension. In addition, a small spurious $m_{\rm stag}$ is found to vanish in the limit of infinite environment dimension. Hence, we propose that this state offers a realization of the QCP or QCPh.

II Symmetric PEPS ansätze

We wish here to consider transitionally invariant fully symmetric PEPS in order to (i) reduce the number of independent variational parameters and (ii) provide a good description of the critical point (or phase) where both SU(2) and lattice symmetries are preserved. For this purpose, we shall use the elegant classification of SU(2)-invariant PEPS tensors on the square lattice Mambrini et al. (2016) according to (i) their virtual degrees of freedom and (ii) how they transform w.r.t the (lattice) point group symmetries (see Fig. 2(a)). For simplicity, we shall a priori restrict ourselves to tensors fully invariant under all operations of the $C_{4v}$ point group (i.e. belonging to the so-called $A_{1}$ IRREP). The tensors are further classified according to their virtual space $V$ given by a direct sum of SU(2) IRREPs or “spins”, i.e. $V=\bigoplus_{\alpha}s_{\alpha}$. We restrict hereafter to bond dimension $D\leq 7$. Among all the possible cases listed in Table 1 we focus on the most interesting ones carrying low virtual spins defined by $V=\frac{1}{2}\oplus 0$ ($D=3$), $V=\frac{1}{2}\oplus 0\oplus 0$ ($D=4$), $V=1\oplus\frac{1}{2}$ ($D=5$), $V=\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ ($D=5$), $V=\frac{1}{2}\oplus\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ ($D=7$), spanned by a small number $\cal D$ of independent tensors, ${\cal D}=2,8,4,10,30$ respectively, given in the Supplementary Materials of Ref. Mambrini et al., 2016 (except for $D=7$ given in the Supplementary Materials of this paper Poilblanc and Mambrini ). Note that a $\pi$-rotation of the spin basis is assumed on the sites of one of the two sublattices of the square lattice. In this basis, a genuine ${\bf q}={\bf q}_{\rm AF}\equiv(\pi,\pi)$ (spontaneous) magnetic order translates into a uniform ${\bf q}=0$ (spontaneous) magnetization. Subsequently, the generator of $SU(2)$ become invariant only up to translations that map the sublattices to themselves (i.e. shifts over two sites).

The iPEPS method combined with full tensor optimization – We shall now focus on the $J_{1}-J_{2}$ spin-$\frac{1}{2}$ Heisenberg model with NN and NNN antiferromagnetic coupling $J_{1}$ and $J_{2}$, respectively, which we have studied at $J_{2}=0$ in the absence of frustration and, for strong frustration, at $J_{2}=0.5$ and $J_{2}=0.55$. Our first goal is to optimize the variational energy within each $\cal D$-dimensional class of SU(2)-invariant PEPS i.e finding the optimum linear superposition of the $\cal D$ independent tensors of each class. Since the number of variational parameters remains small (maximum of ${\cal D}=30$ for $D=7$) we have used a ”brute force” CG optimization as e.g. given in Numerical Recipes Press et al. (2007). However, this requires an efficient iPEPS computation of the variational energy for any set of variational parameters to “feed” the CG routine. This is performed constructing a self-consistent environment around an active $2\times 2$ cluster (see Fig. 2(b)) using an iterative CTMRG algorithm Nishino and Okunishi (1996); Orús and Vidal (2009); Orús (2012) optimized for spatially symmetric tensors. Indeed, we have introduced simple modifications: (i) we use a unique CTM $C$ tensor (side tensor $T$) which is the same for all corners (edges) and (ii) the basic Singular Value Decomposition (SVD) in each CTMRG step to construct the environment is replaced by a (more stable) ED, the CTM being here a symmetric matrix. The largest environment dimension we could handle was $\chi=400$ and $\chi=294$ for $D=5$ and $D=7$, respectively, for which up to 350 or 400 iterations became necessary to converge the environment. Note that the initial $C$ ($T$) tensor is obtained from the $\mathbb{E}$ tensor of Fig. 2(b) by summing over all external $l$ and $u$ ($u$) indices.

Energetics – Variational energies (per site) in each class of tensors are shown in Fig. 3(a) for $J_{2}=0.5$, as a function of the inverse of the environment dimension $\chi$. A rapid comparison between the different classes (for intermediate $\chi$) reveals that, for identical bond dimension $D$, the classes $V=\frac{1}{2}^{\otimes N}\oplus 0$ with $N=1$, $2$ and $3$ (of bond dimensions $D=3$, $5$ and $7$, respectively) give the best results. Hence, hereafter we shall focus on this PEPS family defined in terms of $N$ “colors” of spin-$\frac{1}{2}$. Note that the case $N=1$ was studied previously in Ref. Wang et al. (2013). Tensors are fully optimized up to a maximum bond dimension $\chi_{\rm opt}$, e.g. for $D=5$, $\chi_{\rm opt}=4D^{2}=100$ and, for $D=7$, $\chi_{\rm opt}=2D^{2}=98$. Then, using environment dimensions $\chi>\chi_{\rm opt}$ together with the fixed optimized tensor obtained at $\chi=\chi_{\rm opt}$, one gets true upper bounds of the variational energy. In contrast, for $D=3$, $\chi_{\rm opt}=12D^{2}=108$ already gives the absolute best tensor with enough accuracy. Generically, we found that the energy always decreases with increasing $\chi$ and, at large enough $\chi$, linear fits can be performed in $1/\chi$ to provide $\chi\rightarrow\infty$ extrapolations, also upper bounds of the ($D$-dependent) variational energies. Note that our $D=7$ extrapolation $-0.49502$ lies within only $0.2\%$ of the extrapolated value $-0.4958$ obtained using cluster update finite size $D=9$ PEPS Wang et al. (2016). We have plotted our ($\chi\rightarrow\infty$) results as a function of $1/D$ in Fig. 3(b) showing perfect consistency with the above-mentioned $D=9$ result together with the DMRG extrapolation $-0.4968$ of Ref. Gong et al., 2014. This agreement is remarkable considering the fact that we use only a unique tensor parametrized by a small number of coefficients. Good variational energies have also been found for the simple NN Heisenberg model ($J_{2}=0$) as well as for larger frustration $J_{2}=0.55$ as shown in Appendix A. Our results are summarized in Table 2 and compared to the best estimates, from Quantum Monte Carlo at $J_{2}=0$ Sandvik and Evertz (2010); Sandvik (2010) and from DMRG Gong et al. (2014), VMC Hu et al. (2013) and finite-size PEPS Wang et al. (2016) at $J_{2}=0.5$ and $J_{2}=0.55$. We note however that our variational energies for $J_{2}=0$ and $J_{2}=0.55$ are slightly less accurate as for $J_{2}=0.5$. In fact, we believe $J_{2c}$ is close to 0.5 and we argue below that our (optimized) PEPS is capable of picking up the critical nature of the QCP or QCPh. For $J_{2}=0.55$ translation symmetry breaking is likely to occur spontaneously, which is not captured by our homogeneous ansatz. The ansatz does not either sustain magnetic LR order, that may explain its lower accuracy at $J_{2}=0$.

III Correlation functions

Once the PEPS |\Psi_{0}\big{>}=|\Psi(D,\chi_{\rm opt})\big{>} has been optimized using the largest possible environment dimension $\chi=\chi_{\rm opt}(D)$, various correlation functions can be computed (e.g. along the ${\bf e}_{x}$ horizontal direction), like (i) the spin-spin correlations,

[TABLE]

(ii) the (connected) longitudinal dimer-dimer correlations,

[TABLE]

and (iii) the (connected) transverse dimer-dimer correlations,

[TABLE]

where dimer operators $D_{\bf i}^{x}={\bf S_{i}}\cdot{\bf S}_{{\bf i}+{\bf e}_{x}}$ and $D_{\bf i}^{y}={\bf S_{i}}\cdot{\bf S}_{{\bf i}+{\bf e}_{y}}$ are oriented either along the ${\bf e}_{x}$ (horizontal) or ${\bf e}_{y}$ (vertical) directions, respectively, and the expectation values are taken in the optimized |\Psi_{0}\big{>} PEPS.

The calculations of correlators are accomplished using the set-up shown in Fig. 4(a-c). Appropriate transfer matrices are used so that one can construct arbitrarily long strips. Here the site tensor is fixed to its optimized output using $\chi=\chi_{\rm opt}(D)$ (hereafter we use $\chi_{\rm opt}=49$ for $D=7$) while the environment dimension $\chi>\chi_{\rm opt}(D)$ can be then further increased to reach convergence, which is easily achieved for short distance $r$. A comparison between the results obtained with the two ansätze $V=\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ (a) and $V=\frac{1}{2}\oplus\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ (b) is shown in Fig. 5 for $J_{2}=0.5$. Although a fast decay of the dimer-dimer correlations is seen in both cases, the behavior of the (staggered) spin-spin correlations is qualitatively different : for $D=7$ $|C_{\rm s}(r)|$ seems to approach a finite value while, for $D=5$ (or $D=3$ as well), it steadily decays to zero. This signals the emergence, for $D\geq 7$, of a finite staggered magnetization as defined by $m_{\rm stag}(\chi)=\sqrt{\lim_{r\rightarrow\infty}|C_{\rm s}(r)|}\,.$ We note however that, strictly speaking, for finite $\chi$ the above limit should vanish since the correlations are cut-off above some correlation length $\xi_{\rm s}(\chi)$ (see below). In other words, the strip of Fig. 4(a) is, crudely speaking, similar to a quasi-1D physical strip (ladder) of effective width $L_{\rm eff}(\chi)$ Nishino et al. (1996), which can not sustain long-range magnetic order from Mermin-Wagner theorem (MWT) Mermin and Wagner (1966). However, MWT may not, strictly speaking, apply to a transfer operator as for a true Hamiltonian. In addition, for $D=7$ the SU(2) symmetry is spontaneously broken : small deviations from a perfectly SU(2)-symmetric environment act as a small symmetry-breaking (AF) ”field” and the local spin operator acquires a finite value \big{<}{\bf S_{i}}\big{>}_{0}=\cos{({\bf q}_{\rm AF}\cdot{\bf i})}\,{\bf m}_{\rm stag} oscillating at the antiferromagnetic wave vector ${\bf q}_{\rm AF}$. As shown in Appendix B, $m_{\rm stag}(\chi)$ vanishes in the $\chi\rightarrow\infty$ limit, physically corresponding to the limit of an infinitely wide strip $L_{\rm eff}\rightarrow\infty$. This implies that the infinite 2D system recovers the full SU(2) spin symmetry encoded in the tensor ansatz. We have seen similar behaviors also for $J_{2}=0$ and $J_{2}=0.55$ as well (see Appendix B). Interestingly, the scaling of $m_{\rm stag}$ to zero may depend slightly of the initial CTM of the CTMRG procedure to converge the environment. In contrast, for $D=3$ and $D=5$ the system remains spin isotropic even for finite $\chi$, the spin correlators \big{<}S_{\bf i}^{\alpha}S_{\bf j}^{\alpha}\big{>}_{0} being independent on $\alpha=x,y,z$, as checked explicitly. This signals a qualitative change of behavior when $N\geq 3$ which we identify in the next section.

Diverging correlation lengths – The results described above give some hints that, when $D=7$, the spin-spin correlations become algebraic at long distance. However, for finite bond dimension $\chi$, the strips of Fig. 4(a-c) can be seen as effective 1D systems. Then, finite correlation lengths $\xi_{D}(\chi)$ naturally emerge as the inverse of the gaps of finite-dimensional $D_{\rm eff}^{2}\times D_{\rm eff}^{2}$ transfer matrices, where $D_{\rm eff}=D\chi$ (Fig. 4(a,b)) or $D_{\rm eff}=D^{2}\chi$ (Fig. 4(c)) are the effective dimensions of the associated 1D MPS. Using empirical findings for the correlation length $\xi_{\rm 1D}$ in critical 1D systems Tagliacozzo et al. (2008); Pollmann et al. (2009); Pirvu et al. (2012), $\xi_{\rm 1D}(D)\sim D^{\kappa}$, one then expects that $\xi_{D}(\chi)\sim(D_{\rm eff})^{\kappa}$, $\kappa>0$, which should diverge with $\chi$ as a power law for critical PEPS. Hence, criticality (if any) is restored only in the $\chi\rightarrow\infty$ limit and finite-$\chi$ scaling is necessary to obtain informations on the QCP or QCPh. Note that, when spin rotational symmetry is (artificially) broken at finite $\chi$, it is important to consider the connected spin-spin correlator ${\tilde{C}}_{\rm s}(d)=C_{\rm s}(d)-(m_{\rm stag})^{2}$. From straightforward fits of the long-distance correlations at $J_{2}=0.5$ (see Appendix C) we have extracted the correlation lengths $\xi_{D}(\chi)$ associated to the ${\tilde{C}}_{\rm s}$, $C_{\rm d}^{(\rm T)}$ and $C_{\rm d}^{(\rm L)}$ correlation functions and results are shown in Fig. 6. For $D=3$ or $D=5$ we find a clear saturation of the spin-spin correlation lengths to small values while the dimer-dimer correlations lengths diverge linearly with $\chi$. Such a behavior is typical of bi-partite dimer models Sandvik and Moessner (2006) or of the NN RVB state on the square lattice Albuquerque and Alet (2010); Poilblanc et al. (2012) due to $U(1)$-gauge symmetry. In fact, the $D=3$ PEPS can be viewed as an extended-range RVB state Wang et al. (2013) and the $D=5$ PEPS as an extended-range two-color RVB state. Plotting the dimer correlation lengths in Fig. 6(c,e) as a function of $\chi/D^{2}$ clearly reveals the similarities between $D=5$ and $D=7$. However, in the case of the spin correlations, a sudden qualitative change occurs at $D=7$ for which we find that the spin-spin correlation length no longer saturates but increases linearly with $\chi$, as the dimer correlation lengths do (see Appendices C and D for details). No sign of saturation of the correlation lengths is observed up to the largest available environment dimensions. This suggests that the (optimized) $D=7$ PEPS is critical in the limit $\chi\rightarrow\infty$ or, at least, can very well describe a critical state.

Power-law exponents – 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<d<\xi_{D}$, where $\eta_{\rm s}$ and $\eta_{\rm d}$ defined e.g. in Ref. Lou et al. (2009) are the anomalous dimensions. Note however that this scaling regime can be observed only when $\xi_{D}(\chi)$ has reached a sufficiently large value. To obtain estimates of the exponents $1+\eta_{\rm s}$ and $1+\eta_{\rm d}$ we have plotted spin-spin and (longitudinal) dimer-dimer correlations at $J_{2}=0.5$ in Fig. 7(a,b) using log-log scales. For $D=3$ ($D=5$) the dimer correlation length is very large (is large) for the largest $\chi$ we can reach and, from fits of the data in the range $1<d<100$ ($1<d<20$), one can easily extract the exponent $1+\eta_{\rm d}\simeq 1.25$ ($1+\eta_{\rm d}\simeq 1.5$). For $D=7$, it is difficult to extract accurate exponents since cross-overs to exponential decays occur rapidly around $d\sim\xi_{7}\simeq 6$, for both the spin-spin and dimer-dimer correlations. However, the systematic trend of the data with $\chi$ in Fig. 7(a,b) suggests $\eta_{\rm s}\sim 0.6$ and $\eta_{\rm d}\sim 1.2$.

Discussion and outlook – Above, we have found solid evidence that the $N=3$ ($D=7$) SU(2)-invariant state exhibits slowly (possibly power law) decaying spin-spin and dimer-dimer correlation functions, suggesting a critical behavior, or at least very large correlation lengths. We now argue that the family of SU(2)-symmetric tensors characterized by the virtual space $V=\frac{1}{2}^{\otimes N}\oplus 0$ with $N\geq 3$ “colors” can describe faithfully the QCP or QCPh of the spin-$\frac{1}{2}$ $J_{1}-J_{2}$ Heisenberg model.

First, we observed that spin-spin correlations decay less and less rapidly for increasing $N$ (i.e. $D$) so we expect such correlations to become longer and longer range for increasing $N$. Since the anomalous dimension $\eta_{\rm s}$ (defined from the correlation at intermediate distances) generically decreases with increasing $D$, one can put an upper bound to its infinite-$D$ limit, namely $\eta_{\rm s}<0.6$.

Secondly, it is remarkable that dimer-dimer correlations (and correlation lengths) become very similar for $N=2$ and $N=3$, if compared at the same value of the ratio $\chi/D^{2}$. In fact, we may speculate that, for $N\geq 3$, all correlation lengths diverge as

[TABLE]

where the prefactor $f_{D}$ depends weakly on $D$, the main effect of increasing the bond dimension being to rescale the environment dimension $\chi\rightarrow\chi_{D}=\chi/D^{2}$. We note nevertheless that, although our data are consistent with (6), one cannot rule out that some of the correlation lengths may saturate to a finite, although large, value.

Related $J-Q$ models can be investigated with QMC Lou et al. (2009) and $\eta_{\rm s}\simeq 0.35(2)$ and $\eta_{\rm d}\simeq 0.20(2)$ have been obtained (for the $J-Q_{2}$ model), which seem to deviate substantially from our estimates above. However, our estimation of $\eta_{\rm s}$ seems consistent with the VMC result Morita et al. (2015) $\eta_{\rm s}\sim 0.5$ obtained for the $J_{1}-J_{2}$ Heisenberg model at $J_{2}=0.5$.

Note that the power-law exponent $1+\eta_{\rm d}$, extracted from the correlations at intermediate distances $d<\xi_{D}(\chi)$, seems to increase significantly with $D$. The predicted large value of the $D\rightarrow\infty$ dimer anomalous dimension might indicate that dimer correlations at the QCP or within the QCPh are significantly suppressed compared to $J-Q$ models.

Acknowledgements.

This project is supported by the TNSTRONG ANR grant (French Research Council). This work was granted access to the HPC resources of CALMIP supercomputing center under the allocations 2016-P1231 and 2017-P1231. DP thanks Nicolas Renon (CALMIP) and Cyril Mazauric (ATOS) for assistance. D.P. acknowledges illuminating discussions with Philippe Corboz, as well as helpful advices to implement the CTM algorithm. M.M. and D.P. thank Roman Orus for insightful comments. D.P. also acknowledges inspiring conversations with Fabien Alet, Federico Becca, Ignacio Cirac, Shenghan Jiang, Naoki Kawashima, Frédéric Mila, David Perez-Garcia, Frank Pollmann, Pierre Pujol, Ying Ran, Anders Sandvik, Norbert Schuch, Frank Verstraete and Ling Wang.

Appendix A Scaling of the $D=7$ variational energy vs inverse environment dimension

We report in Fig. 8(a-c) the variational energies of the $V=\frac{1}{2}\oplus\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ PEPS ansatz for the $J_{1}-J_{2}$ model at $J_{2}=0$ (unfrustrated case), $J_{2}=0.5$ and $J_{2}=0.55$. The parameters of the PEPS are optimized with an environment dimension $\chi_{\rm opt}=D^{2}=49$, independently for each value of $J_{2}$. For $J_{2}=0.5$ we also carried out the optimization with $\chi_{\rm opt}=2D^{2}=98$, providing a slightly better energy. The environment dimension $\chi>\chi_{\rm opt}$ is then increased, keeping the PEPS tensor fixed, and the energy is extrapolated linearly with $1/\chi$. At $J_{2}=0.5$, an excellent agreement is found with extrapolation from $D=9$ PEPS cluster update Wang et al. (2016). For $J_{2}=0$ and $J_{2}=0.55$ a less good agreement is found with QMC Sandvik and Evertz (2010); Sandvik (2010) and $D=9$ PEPS cluster update Wang et al. (2016), respectively (see text for explanation).

Appendix B Scaling of the $D=7$ staggered magnetization vs inverse environment dimension

We report in Fig. 9(a-c) the spurious staggered magnetization of the $V=\frac{1}{2}\oplus\frac{1}{2}\oplus\frac{1}{2}\oplus 0$ PEPS ansatz for the $J_{1}-J_{2}$ model at $J_{2}=0$ (unfrustrated case), $J_{2}=0.5$ and $J_{2}=0.55$ (optimized using $\chi_{\rm opt}=D^{2}=49$). The procedure is the same as in Appendix A and the data are plotted vs $\chi$. For all $J_{2}$ values, the scaling (algebraic fits) is consistent with vanishing $m_{\rm stag}$ when $\chi\rightarrow\infty$. Full SU(2) invariance is recovered in this case.

Appendix C Extracting the correlation lengths $\xi_{D}(\chi)$ from the long distance correlations

In order to extract the correlation lengths associated to the various correlation functions $C_{\lambda}(d)$ ($\lambda={\rm S,D}$) defined in the paper in Eqs. (1), (2) and (3), we have computed the long-distance correlations using the transfer matrix methods sketched in Fig. 2. Due to a finite gap in the relevant transfer matrices for all finite dimensions $D$ and $\chi$, one expects an exponential decay of all correlations,

[TABLE]

at sufficiently large distance $d$ (typically $d>\xi_{D}(\chi)$). Let us summarize the procedure : First, the local tensors for $D=3$, $5$ and $7$ are obtained by a full CG optimization (for $J_{2}=0.5$) using a given environment dimension $\chi_{\rm opt}=108$, $100$ and $49$, respectively. The correlations in these fixed PEPS are then computed for increasing values of the environment dimension $\chi$ in two steps : (i) For every choice of $\chi\geq\chi_{\rm opt}$, the new converged CTM $C$ and edge tensor $T$ are computed (by the iterative renormalization scheme) and, finally, (ii) used to compute the correlation functions in the setup shown in Fig. 2(a-c). Results are displayed using semi-logarithmic scales in Figs. 10(a), 11(a) and 12(a). By fitting the asymptotic linear behaviors of the data according to $\ln{C_{\lambda}(d)}=-(1/\xi)d+c_{0}$, one straightforwardly gets the correlation lengths $\xi$ from the slopes $-1/\xi$.

The scaling of the correlation lengths $\xi_{D}$ with $\chi$ are shown in Figs. 10(b), 11(b,c) and 12(b,c). For $D=3$ and $D=5$, one observes a clear saturation of the spin correlation lengths $\xi_{3}$ and $\xi_{5}$ to rather small values (less than 2 lattice spacings) while the dimer correlation length scales linearly with $\chi$ suggesting that $\xi_{D}\rightarrow\infty$ in the limit $\chi\rightarrow\infty$, for which the calculation becomes exact. Note that the (extrapolated) spin correlation length increases with $D$ while the divergence of the dimer correlation length becomes weaker. For $D=7$, one has to consider the connected part of the spin-spin correlation, subtracting off the contribution from the spurious staggered spin density background. The spin correlation length no longer saturates but rather increases linearly with the environment dimension $\chi$. This strongly suggests that $\xi_{7}$ diverges in the limit $\chi\rightarrow\infty$, that is consistent with a power-law decay of the correlation function. We believe our numerical results also support the divergence of both dimer-dimer correlation lengths, as well. Note however that, although the transverse and longitudinal dimer-dimer correlation lengths seem to match for $D=3$ and $D=5$, they deviate substantially for $D=7$, which may be related to the non-vanishing of the spin-spin correlation in that case.

Appendix D Comparison between correlation functions in the $D=7$ PEPS

In principle, correlation lengths can also be extracted directly from the low-energy eigenvalues of the zero dimensional transfer matrix of the one-dimensional tensor network structures arising in Fig. 4. It would be the same transfer matrix for spin-spin and (longitudinal) dimer-dimer correlation function, but the difference would be how the corresponding virtual eigenvectors of these eigenvalues transform under the symmetry. In a perfectly SU(2)-symmetric state giving rise to a SU(2)-symmetric environment (as it occurs for $D=3$ and $D=5$), different selection rules for the singlet (dimer) and the triplet (spin) channels lead to separate blocks of the transfer matrix and, hence, to different correlation lengths, in agreement with our findings. However, for $D=7$ spontaneous SU(2) symmetry breaking occurs and the environment acquires some (staggered) magnetization $\bf m_{\rm stag}$. We believe spin-rotational invariance ($U(1)$ symmetry) is still preserved around the direction of the staggered magnetization. The latter can be pointing in any (arbitrary) direction in the $(x,z)$ plane making difficult the symmetry analysis of the zero dimensional transfer matrix arising in Fig. 4. Analysis of the correlation functions given e.g. by Eqs. 1, 2 or 3 is more straightforward.

At this point, it is not clear whether the long distance spin correlation described in the text is an artifact of the symmetry breaking that i) may lead to a mixture of (diverging) singlet and (short-range) triplet correlations or ii) may lead to ”Goldstone critical behavior” of the transverse spin correlation function. We give arguments below that none of the above applies and argue that the critical behavior of the spin correlation function is an intrinsic feature of the $D=7$ PEPS spin liquid.

For this purpose, we decompose the local spin operator into its longitudinal and transverse spin components,

[TABLE]

where $\bf n$ is a unit vector along ${\bf m}_{\rm stag}$, $S_{\bf i}^{\parallel}={\bf S}_{\bf i}\cdot{\bf n}$ and ${\bf S}_{\bf i}^{\perp}={\bf S}_{\bf i}-({\bf S}_{\bf i}\cdot{\bf n})\,{\bf n}$. The spin correlation function can be then split into its longitudinal and transverse components as $C_{\rm s}(d)=C_{\rm s}^{\parallel}(d)+C_{\rm s}^{\perp}(d)$ with,

[TABLE]

For a true singlet wave function (for which ${\bf m}_{\rm stag}={\bf 0}$), whatever the choice of the vector $\bf n$, on gets $C_{\rm s}^{\perp}(d)=2C_{\rm s}^{\parallel}(d)$. As shown in Fig. 13(a) this is also true for the $D=7$ PEPS, at short distance only (in semi-log scale the two curves are just shifted by $\ln{2}$). At longer distance, however, the longitudinal and transverse spin correlations show different exponential decays. As shown in Fig. 13(b) the correlation length of the longitudinal correlations is much shorter than the one of the transverse correlations. However, both seem to diverge with increasing $\chi$, suggesting that both correlators are critical, possibly power-law, in the $\chi\rightarrow\infty$ limit. This is different from a ”Golstone mechanism” for which the longitudinal correlations remain short-range. Finally, we compare the two spin correlation lengths to the (longitudinal) dimer correlation length. Fig. 13(b) shows that none of the three (diverging) correlation length match, suggesting that the (supposedly) critical behavior of the spin-spin correlation is not induced by the critical behavior of the dimer correlation and is an intrinsic property of the $D=7$ PEPS.