Tensor network study of the spin-1/2 square-lattice $J_{1}$ - $J_{2}$ - $J_{3}$ model: incommensurate spiral order, mixed valence-bond solids, and multicritical points

Wen-Yuan Liu Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA Didier Poilblanc Laboratoire de Physique Théorique, C.N.R.S. and Université de Toulouse, 31062 Toulouse, France Shou-Shu Gong School of Physical Sciences, Great Bay University, Dongguan 523000, China, and
Great Bay Institute for Advanced Study, Dongguan 523000, China Wei-Qiang Chen Shenzhen Institute for Quantum Science and Engineering and Department of Physics, Southern University of Science and Technology, Shenzhen 518055, China Shenzhen Key Laboratory of Advanced Quantum Functional Materials and Devices, Southern University of Science and Technology, Shenzhen 518055, China Zheng-Cheng Gu Department of Physics, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong, China

(January 27, 2025; January 27, 2025)

Abstract

We use the finite projected entangled pair state (PEPS) method to investigate the global phase diagram of the spin-1/2 square-lattice $J_{1}$ - $J_{2}$ - $J_{3}$ antiferromagnetic (AFM) Heisenberg model. The ground-state phase diagram is established with a rich variety of phases: Néel AFM, gapless quantum spin liquid, valence-bond solid (VBS), stripe AFM, and incommensurate spiral phases. The nature of the VBS region is revealed, which contains a plaquette VBS and a mixed columnar-plaquette VBS, with the emergence of short-range incommensurate spin correlation in some region. The long-range incommensurate magnetic phase is also explicitly characterized as a planar spiral with incommensurate spatial periodicities. Most interestingly, there exists several multicritical points connecting different phases. These findings elucidate the true nature of the long-standing square-lattice $J_{1}$ - $J_{2}$ - $J_{3}$ antiferromagnet at zero-temperature. Our results also pave the way to accurately simulate complex two-dimensional quantum spin systems that may host nonuniform features by means of the finite PEPS.

I Introduction

The spin-1/2 $J_{1}$ - $J_{2}$ - $J_{3}$ antiferromagnetic (AFM) Heisenberg model on the square lattice is one of the paradigmatic prototypes to study frustrated quantum magnets. This model has attracted a lot of interests after the application of Anderson’s resonating valence-bond (RVB) theory to high-temperature superconductivity Anderson (1987). The Hamiltonian of this system is described as follow:

H=J_{1}\sum_{\langle i,j\rangle}\mathbf{S_{i}}\cdot\mathbf{S_{j}}+J_{2}\sum_{\langle\langle i,j\rangle\rangle}\mathbf{S_{i}}\cdot\mathbf{S_{j}}+J_{3}\sum_{\langle\langle\langle i,j\rangle\rangle\rangle}\mathbf{S_{i}}\cdot\mathbf{S_{j}},\quad

(1)

where $J_{1}$ , $J_{2}$ , $J_{3}$ denote the first-, second- and third-nearest neighbour couplings respectively, and the summations run over all corresponding spin pairs. For the positive $J_{i}$ couplings, the three kinds of interactions compete with each other and thus leads to intractable difficulties for analytic and numerical studies. The classical phase diagram obtained by spin-wave theory contains four phases: (i) a Néel AFM phase ordered at wave vector ${\bf k_{0}}=(\pi,\pi)$ , (ii) a stripe AFM phase ordered at wave vector ${\bf k_{x}}=(\pi,0)$ or ${\bf k_{y}}=(0,\pi)$ , (iii) a spiral phase at ${\bf Q_{x}}=(\pm q,\pi)$ or ${\bf Q_{y}}=(\pi,\pm q)$ with $\cos{q}=(2J_{2}-J_{1})/4J_{3}$ , and (iv) another spiral phase at ${\bf Q}=(\pm q,\pm q)$ with $\cos{q}=-J_{1}/(2J_{2}+4J_{3})$ Gelfand et al. (1989); Moreo et al. (1990), as shown in Fig. 1(a). The Néel AFM and the spiral $(\pm q,\pm q)$ phases are separated by a classical critical line $(J_{2}+2J_{3})=J_{1}/2$ . However, the quantum phase diagram with spin- $1/2$ is not fully understood.

For spin- $1/2$ , it has been established that at small $J_{2}$ and $J_{3}$ couplings, the model possesses a Néel AFM order Chandra and Douçot (1988); Figueirido et al. (1990); Moreo et al. (1990); Chubukov (1991). At larger $J_{2}$ and $J_{3}$ , the corresponding stripe AFM and spiral orders will develop Read and Sachdev (1991). The controversy exists in the intermediate region of $J_{2}$ and $J_{3}$ , especially along the classical critical line. Some theories strongly suggest that the combined effect of enhanced quantum fluctuations and frustration could destroy the long-range orders and stabilize paramagnetic states along the critical line, including the spin-wave theories Chandra and Douçot (1988); Moreo et al. (1990); Chubukov (1991), renormalization group (RG) analysis of the non-linear $\sigma$ model Ioffe and Larkin (1988), series expansions Gelfand et al. (1989), and momentum-shell RG calculation Einarsson and Johannesson (1991), etc. There are also other theories supporting the existence of paramagnetic states but for the $J_{2}$ and $J_{3}$ couplings shifted to larger values with respect to the classical critical line by quantum fluctuations Read and Sachdev (1991); Ferrer (1993).

From the combined analyses of different studies, it seems most likely there exists an intermediate nonmagnetic region. Nevertheless, the nature of the nonmagnetic region is far from clear. While the large- $N$ and series expansion results predict that this intermediate state is spontaneously dimerized Read and Sachdev (1991); Gelfand et al. (1989), the spin-wave theory suggests a spin liquid state Chandra and Douçot (1988); Moreo et al. (1990); Chubukov (1991). For the $J_{1}$ - $J_{3}$ model with $J_{2}=0$ , a Monte Carlo study of the classical limit supplemented by analytical arguments on the role of quantum fluctuations supports the emergence of a valence-bond solid (VBS) or a $Z_{2}$ spin liquid between the Néel AFM and spiral phases Capriotti and Sachdev (2004). While exact diagonalization calculations on small clusters suggest a VBS state with incommensurate short-range spin correlation Leung and Lam (1996), the early density matrix renormalization group (DMRG) results may support a gapped spin liquid Capriotti et al. (2004). For $J_{2}\neq 0$ , a short-range valence bond study finds a plaquette VBS state along the line with $J_{2}+J_{3}=J_{1}/2$ (where the description in terms of nearest-neighbor singlet coverings is excellent) Mambrini et al. (2006). Later, a mixed columnar-plaquette VBS state was also proposed Ralko et al. (2009), which possesses a long-range plaquette order but breaks the isotropy between the $x$ and $y$ directions, supporting the findings in a previous study du Croo de Jongh et al. (2000). Furthermore, exact diagonalizations also show that quantum fluctuations could lead to new quantum phases Sindzingre et al. (2010).

Refer to caption — Figure 1: (a) Classical phase diagram for the ground state of the $J_{1}$ - $J_{2}$ - $J_{3}$ model: the AFM phase ordered at $(\pi,\pi)$ , the stripe phase ordered at $(0,\pi)$ or $(\pi,0)$ , a spiral phase ordered at $(\pm q,\pm q)$ , and another spiral phase ordered at $(\pi,\pm q)$ or $(\pm q,\pi)$ . (b) Quantum phase diagram for the ground state of the spin-1/2 square-lattice $J_{1}$ - $J_{2}$ - $J_{3}$ AFM model: the AFM phase ordered at $(\pi,\pi)$ , the stripe phase ordered at $(0,\pi)$ or $(\pi,0)$ , the spiral phase ordered at $(\pm q,\pm q)$ , the gapless QSL phase (red) and the VBS phase. Tricritical points (end points) are shown by red (blue) dots. The AFM-spiral, AFM-stripe, VBS-spiral and VBS-stripe transitions are first-order, and other phase transitions are continuous. The phase transition lines are determined via finite size extrapolations. In the VBS region, below the violet dashed line is a plaquette VBS and above is a mixed columnar-plaquette VBS. The region in the right corner might contain a quadruple point (the white dot) connecting the VBS, QSL, AFM and stripe phases, or two multicritical points shown in (c) (assuming the gapless QSL does not touch the stripe phase). Grey dashed-dotted lines marking $J_{2}=0$ and $J_{3}=0$ are shown in (b) and (c).

The special case of the $J_{1}$ - $J_{2}$ - $J_{3}$ model with $J_{3}=0$ , i.e., the $J_{1}$ - $J_{2}$ model, is of great interest in the field of quantum magnetism, which has been intensively studied for more than three decades. The $J_{1}$ - $J_{2}$ model exhibits a Néel AFM phase in the region $0\leq J_{2}/J_{1}\lesssim 0.45$ and a stripe AFM phase for $J_{2}/J_{1}\geq 0.61$ . The intermediate nonmagnetic region has been investigated by many different methods Chandra and Douçot (1988); Dagotto and Moreo (1989); Figueirido et al. (1990); Sachdev and Bhatt (1990); Poilblanc et al. (1991); Chubukov and Jolicoeur (1991); Schulz and Ziman (1992); Ivanov and Ivanov (1992); Einarsson and Schulz (1995); H. J. Schulz and Poilblanc (1996); Zhitomirsky and Ueda (1996); Singh et al. (1999); Capriotti and Sorella (2000); Capriotti et al. (2001); Zhang et al. (2003); Takano et al. (2003); Sirker et al. (2006); Schmalfuß et al. (2006); Mambrini et al. (2006); Darradi et al. (2008); Arlego and Brenig (2008); Isaev et al. (2009); Murg et al. (2009); Beach (2009); Richter and Schulenburg (2010); Yu and Kao (2012); Jiang et al. (2012); Mezzacapo (2012); Wang et al. (2013); Hu et al. (2013); Doretto (2014); Qi and Gu (2014); Gong et al. (2014); Chou and Chen (2014); Morita et al. (2015, 2015); Richter et al. (2015); Wang et al. (2016); Poilblanc and Mambrini (2017); Wang and Sandvik (2018); Haghshenas and Sheng (2018); Liu et al. (2018); Poilblanc et al. (2019); Hasik et al. (2021); Ferrari and Becca (2020); Liu et al. (2022a); Nomura and Imada (2021), and the existence of two successive gapless quantum spin liquid (QSL) and VBS phases lying between the Néel and stripe phase is being supported by more recent works Wang and Sandvik (2018); Ferrari and Becca (2020); Liu et al. (2022a); Nomura and Imada (2021). In particular, the finite PEPS method provides very solid large-size results and makes it possible to extract critical exponents to understand the physical nature of the AFM-QSL and QSL-VBS phase transitions Liu et al. (2022a), which indicates an intrinsic relation between the gapless QSL and deconfined quantum critical point (DQCP) Senthil et al. (2004a).

Very recently, a tensor network study has also revealed the exotic properties of the $J_{1}$ - $J_{2}$ - $J_{3}$ model Liu et al. (2022b). Specifically, the highly accurate tensor network results from finite and infinite PEPS simulations not only identify the extended QSL and VBS phases at finite $J_{3}$ , but also demonstrate the intrinsic relation between the gapless QSL and DQCP, offering a fantastic scenario to understand the two exotic quantum phenomena. The existence of QSL and VBS phases in the pure $J_{1}$ - $J_{3}$ model, i.e., $J_{2}=0$ , is also supported by recent DMRG calculations Wu et al. (2022).

However, because of the highly complex competition of $J_{1}$ , $J_{2}$ and $J_{3}$ interactions, the global phase diagram of the $J_{1}$ - $J_{2}$ - $J_{3}$ model remains enigmatic, particularly in the regions that may involve nonuniform physical properties like the incommensurate spin correlations. Thanks to the advancement of tensor network method, we can now study this model in a precise way. One approach to study this model is based on the widely used infinite PEPS (iPEPS) ansatz Liu et al. (2022b). Although it has a presumed unit cell, the nonuniform properties can still be investigated by systematically enlarging the size of unit cell. There have been remarkable progresses in correlated electron models along this direction Corboz et al. (2011, 2014); Zheng et al. (2017). In this paper, we shed new light on the ground-state phase diagram of the $J_{1}$ - $J_{2}$ - $J_{3}$ model by using an alternative approach, namely, the finite PEPS ansatz, which is well-suited for describing nonuniform features with site-independent tensors. The finite PEPS works very well in the scheme combined with variational Monte Carlo, and the full details of the algorithm can be found in the preceding publications Liu et al. (2017, 2021, 2022a, 2022b). The obtained phase diagram is presented in Fig. 1(b).

Essentially, this paper is an extension of our previous work on the $J_{1}$ - $J_{2}$ - $J_{3}$ model Liu et al. (2022b), but emphasizing different aspects. In the previous work, we discovered a large region of gapless QSL between the Néel AFM and VBS phases. By tuning coupling constants we found that the extension of the gapless QSL decreases and eventually merges into a line of direct continuous transition between the Néel AFM and VBS phases, showing a direct connection between gapless QSL and DQCP. With further analyses of the AFM-VBS, AFM-QSL and QSL-VBS transitions, we explicitly demonstrated a new scenario to understand gapless QSL and DQCP, suggesting they are described by a unified quantum field theory. The gapless QSL and DQCP physics actually is only a part of the $J_{1}$ - $J_{2}$ - $J_{3}$ model, and here we turn to reveal its global ground state phase diagram. Specifically, we comprehensively investigate its nonuniform features including incommensurate short-range and long-range spiral properties. We also find a novel mixed valence-bond solid phase as well as several multicritical points.

Rest of the paper is organized as follows. In Sec. II, we first present the results obtained at $J_{2}=0$ . By computing spin and dimer order parameters, as well as ground state energies, we provide a complete phase diagram of the $J_{1}$ - $J_{3}$ model, and focus on the emergence of incommensurate spiral spin correlations upon increasing $J_{3}$ . Then, we consider an intermediate $J_{2}\neq 0$ term to further analyze the properties of the VBS phase. Next we consider a large negative $J_{2}$ , which can suppress the intermediate VBS phase and eventually lead to a direct first-order AFM-spiral phase transition. Finally, we consider the $J_{3}<0$ case to obtain the complete phase diagram, where the gapless QSL and VBS phases eventually disappear, and instead a direct AFM-stripe transition occurs. In Sec. III, we discuss the existing results and obtain a refined description of the phase diagram of the $J_{1}$ - $J_{2}$ - $J_{3}$ model. Finally, three appendices are provided, including the convergence with bond dimension, the ED results of a $4\times 4$ cluster, and a refined analysis of the spin structure factor in the mixed VBS phase.

II results

We use the recently developed finite PEPS method to perform all the calculations Liu et al. (2017, 2021). The finite PEPS method works very well in the scheme of variational Monte Carlo approach, where the summation of physical degrees of freedom is replaced by Monte Carlo sampling. Thus, one only needs to deal with the single-layer tensor networks at a computational cost scaling ${\cal{O}}(D^{6})$ , where $D$ is the tensor bond dimension. Such an approach has been successfully applied to solve very challenging quantum many-body problems Liu et al. (2022a, b, c, 2023). In Ref. Liu et al. (2022b), we have demonstrated that the simulations with $D=8$ can obtain well converged results for the $J_{1}$ - $J_{2}$ - $J_{3}$ model up to the $20\times 28$ system size. Here we have checked the $D$ -convergence on the open $16\times 16$ cluster at $(J_{2},J_{3})=(-1,0.8)$ and find that $D=8$ can also provide converged energy and order parameters, as shown in Appendix A. Therefore, we use $D=8$ for all the calculations and we set $J_{1}=1$ as the energy unit unless otherwise specified.

II.1 $J_{2}=0$

II.1.1 Magnetic order parameters

Through detailed computations, we find four phases by varying $J_{3}$ : (I) the Néel AFM phase ordered at ${\bf k_{0}}=(\pi,\pi)$ , (II) the gapless QSL phase, (III) the VBS phase, and (IV) the long-range spiral ordered phase. These four phases have been displayed in Fig. 2 along the $J_{3}$ axis. Note that, when $J_{3}$ is sufficiently large, phase (IV) evolves continuously to approach the commensurate spiral state ordered at ${\bf k_{1}}=(\pm\pi/2,\pm\pi/2)$ .

The gapless spin liquid phase and associated AFM-QSL and QSL-VBS phase transitions have been studied in a previous finite PEPS simulation work Liu et al. (2022b). Now we focus on the VBS-spiral transition. We define two local order parameters to detect long-range spin ordered phases:

	$\displaystyle\langle S_{z}^{\rm abs}\rangle=\frac{1}{L^{2}}\sum_{i_{x},i_{y}}\|\langle S^{z}_{i_{x},i_{y}}\rangle\|,$		(2)
	$\displaystyle\langle S_{z}^{\rm AFM}\rangle=\frac{1}{L^{2}}\sum_{i_{x},i_{y}}(-1)^{i_{x}+i_{y}}\langle S^{z}_{i_{x},i_{y}}\rangle.$		(3)

Since $\langle S_{z}^{\rm abs}\rangle$ is the average of the absolute value of the onsite $\langle S^{z}_{i_{x},i_{y}}\rangle$ , it can distinguish magnetic phases from nonmagnetic ones. Note $\langle S_{z}^{\rm AFM}\rangle$ is the standard Néel AFM order parameter, which can distinguish the Néel AFM phase from other phases including magnetic ones. These two order parameters should be zero theoretically on finite-size systems because of the SU(2) symmetry, but it is expected that for magnetic phases the corresponding ground state may break SU(2) symmetry and show nonzero values for $\langle S_{z}^{\rm abs}\rangle$ or $\langle S^{z}_{i_{x},i_{y}}\rangle$ if the system size becomes large enough.

In Fig. 2, we show the variation of $\langle S_{z}^{\rm abs}\rangle$ and $\langle S_{z}^{\rm AFM}\rangle$ with increasing $J_{3}$ on the $16\times 16$ system size. Indeed, in the Néel AFM phase, $\langle S_{z}^{\rm abs}\rangle$ and $\langle S_{z}^{\rm AFM}\rangle$ are equal and finite. After entering into the nonmagnetic QSL and VBS phases, SU(2) symmetry is almost restored, signaled by the much smaller values of $\langle S_{z}^{\rm abs}\rangle$ and $\langle S_{z}^{\rm AFM}\rangle$ close to zero. When approaching the spiral phase, SU(2) symmetry is broken again and $\langle S_{z}^{\rm abs}\rangle$ enhances sharply, indicating a first-order phase transition to a magnetic phase at $J_{3}\simeq 0.75$ . In addition, we observe that $\langle S_{z}^{\rm AFM}\rangle$ still vanishes, which actually is the result of the incommensurate spin pattern.

II.1.2 Nature of the VBS phase

We also check the dimer order parameters, defined as Liu et al. (2022a)

\langle D_{\alpha}\rangle=\frac{1}{N_{b}}\sum_{{\bf i}}(-1)^{i_{\alpha}}\langle B^{\alpha}_{\bf i}\rangle,

(4)

where $B^{\alpha}_{\bf i}={\bf S}_{{\bf i}}\cdot{\bf S}_{{\bf i}+{\rm e_{\alpha}}}$ is the bond operator between site ${\bf i}$ and site ${\bf i}+{\rm e_{\alpha}}$ along $\alpha$ direction ( $\alpha=x$ or $y$ ). $N_{b}=L(L-1)$ is the total number of counted bonds. As seen in Fig. 2 of the results on a finite $16\times 16$ cluster, the VBS order parameters $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ also a show sudden drop from the VBS phase (III) to the spiral phase (IV), consistent with a first-order phase transition scenario.

Furthermore, we observe that $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ are equal for $J_{3}<0.55$ and start to deviate for $J_{3}>0.55$ , suggesting a transition from the plaquette VBS to the mixed-VBS phase. We believe that the difference between $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ might be induced by short-range spiral correlations (SRSC) (shown in Fig. 8 and discussed later). By analyzing spin correlations, we find that the SRSC appears already at $J_{3}\simeq 0.45$ . This observation suggests that the stronger SRSC at $J_{3}\simeq 0.55$ may induce different VBS order parameters $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ , which is a characteristic feature of the mixed-VBS state.

To examine the size dependence of the dimer order parameters, we present $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ with respect to $1/L$ in Fig. 3. One can find clearly that at $J_{3}=0.8$ the dimer orders tend to vanish in the thermodynamic limit, in contrast to other $J_{3}$ cases inside the estimated VBS phase region $0.37<J_{3}\lesssim 0.75$ . For $J_{3}=0.45$ and $0.5$ , through finite-size scaling up to $20\times 20$ sites, the extrapolated $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ values in the thermodynamic limit are identical, indicating a plaquette VBS order. For $J_{3}=0.6$ and $0.7$ , the extrapolated values of $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ are nonzero but different, suggesting a mixed columnar-plaquette VBS order. More details are provided in Appendix B. The transition point between the plaquette and mixed VBS phase is estimated to be $J_{3}\simeq 0.55$ according to the results of the $16\times 16$ size, which is shown as the violet dotted line in Fig. 2.

We note that the exact diagonalization of the system with twisted boundary conditions on the size up to $32$ sites, marks the region $0.5\lesssim J_{3}\lesssim 0.57$ as a mixed columnar-plaquette VBS phase, which is followed by a SRSC state before entering the long-range spiral ordered phase Sindzingre et al. (2010). This picture is roughly consistent with our results, but we believe that the SRSC state found by exact diagonalization is actually part of the VBS phase. In addition, a short-range valence bond study suggests a plaquette VBS state along the line $J_{2}+J_{3}=J_{1}/2$ Mambrini et al. (2006). In combination with our results, the VBS phase is most likely comprised of a plaquette and a mixed columnar-plaquette phase, and the transition between them is continuous. A similar scenario has been reported in a quantum dimer model on the square lattice, where a mixed columnar-plaquette VBS phase continuously intervenes between a columnar and a plaquette VBS phase Ralko et al. (2008). Qualitatively, we know that the competition between $J_{1}$ and $J_{3}$ couplings tends to induce a spiral order. However, when $J_{3}$ is too small to lead to the SRSC, the isotropy between the $x-$ and $y-$ direction survives, consistent with the existence of a plaquette VBS phase. When $J_{3}$ gets larger but still not large enough to stabilize a spiral long-range order, rotation symmetry breaks down and gives rise to a mixed columnar-plaquette VBS phase.

II.1.3 Energy curve versus $J_{3}$

To have a more comprehensive understanding of the various phases, we also consider the $J_{3}$ dependence of ground-state energy. In Fig. 4, we plot the energy per site at different $J_{3}$ based on the $16\times 16$ system size. The filled circles (red and blue) denote the energy per site obtained from all the sites, and the empty ones (red and blue) are the energy per site obtained from the central $(L-8)\times(L-8)$ sites. One can see that the energy of the $16\times 16$ system has a maximum at $J_{3}\approx 0.52$ , which interestingly is accompanied with the sign change of the third-nearest neighbour terms $\langle{\bf S_{i}}\cdot{\bf S_{j}}\rangle$ , i.e. the third-nearest neighbour $\langle{\bf S_{i}}\cdot{\bf S_{j}}\rangle$ is positive for $J_{3}\lesssim 0.52$ but negative for $J_{3}\gtrsim 0.52$ . To locate the first-order transition point from the VBS to spiral phase from the ground-state energy curve, we first use the ground state at $J_{3}=0.65$ (with initialization by the simple update Jiang et al. (2008)) as an initialization for further optimization to obtain the ground state at $J_{3}=0.68$ . Next, we use the state at $J_{3}=0.68$ to get the optimized state at $J_{3}=0.7$ , and further on to get the optimized states of $J_{3}=0.72$ , $0.75$ and $0.78$ sequentially. Such a process can be viewed as an adiabatic evolution process. On the other hand, we can also perform the reverse process starting from the ground state at $J_{3}=0.8$ and move backwards sequentially to get the optimized states at $J_{3}=0.78$ , $0.755$ , $0.75$ and $0.72$ . Except in rare special cases, first-order transitions generically show hysteresis under such an evolution, i.e. one does not adiabatically evolve into another phase. The transition point can be obtained from the crossing of the two energy evolution curves. In Fig. 4, one can see that the energy curves of the $L\times L$ ( $L=16$ ) system obtained from the two different paths (the red and blue solid circles) show a crossing around $J_{3}\simeq 0.75$ . Based on the energy curves of the central $(L-8)\times(L-8)$ sites which can provide a better estimation of the energy in the thermodynamic limit, the crossing is clearer. Both analyses indicate that the VBS-spiral transition occurs at $J_{3}\simeq 0.75$ , in good agreement with the behaviour of the spin order parameter $\langle S^{\rm abs}_{z}\rangle$ which shows an apparent discontinuity at the same $J_{3}$ value, as shown in Fig. 2. We note that the VBS-spiral transition point can even be estimated via the second derivative of the energy with respect to $J_{3}$ on a periodic $4\times 4$ cluster, which provides the identical result as shown in Appendix C.

II.1.4 Incommensurate spiral order

Now we examine how the spiral order evolves with increasing $J_{3}$ , focusing on the $J_{1}$ - $J_{3}$ model, by computing the spin structure factor defined as

S({\bf k})=\frac{1}{L^{2}}\sum_{\bf{ij}}\langle{\bf S}_{{\bf i}}\cdot{\bf S}_{{\bf j}}\rangle{e}^{i{\bf k}\cdot({\bf i}-{\bf j})}.

(5)

We compute all pairs of spin correlators $\langle{\bf S}_{{\bf i}}\cdot{\bf S}_{{\bf j}}\rangle$ to get the structure factor. In the classical $J_{1}$ - $J_{3}$ model, for $J_{3}\leq 0.25$ , the ground state is a conventional Néel AFM state with magnetic wave vector ${\bf k_{0}}=(\pi,\pi)$ . For $J_{3}/J_{1}>0.25$ , the ground state has a planar incommensurate magnetic order at a wave vector ${\bf Q}=(\pm q,\pm q)$ . With increasing $J_{3}/J_{1}$ , the spiral order is incommensurate except at $J_{3}/J_{1}=0.5$ and $J_{3}/J_{1}\rightarrow\infty$ , and the wave vector ${\bf Q}=(\pm q,\pm q)$ will gradually move from $q=\pi$ to $q=\pi/2$ according to the formula $q={\rm cos}^{-1}(-0.25J_{1}/J_{3})$ for $J_{3}/J_{1}>0.25$ Locher (1990); Capriotti and Sachdev (2004).

In Fig. 5, we show the spin structure factor of the $16\times 16$ system size for different $J_{3}$ at $J_{2}=0$ . At $J_{3}=0.2$ (in the Néel AFM phase), $0.35$ (in the QSL phase), and $0.42$ (in the VBS phase), the peaks of spin structure factor are located at ${\bf k_{0}}=(\pi,\pi)$ . Further increasing $J_{3}$ , the peak will gradually move from ${\bf k_{0}}=(\pi,\pi)$ to ${\bf Q}=(\pm q,\pm q)$ . If $J_{3}/J_{1}\rightarrow\infty$ (set $J_{1}=0$ ), the wave vector will approach to ${\bf Q}=(\pm\pi/2,\pm\pi/2)$ , as shown in Fig. 5(h). In Fig. 6, we present the classical and quantum results of the change of wave vector for a clear comparison. We would like to mention that a careful inspection of the data at $J_{3}=0.7$ , as shown in Fig. 5(f), reveals a detectable difference in the magnitudes of the maxima (see Appendix D), which is consistent with the mixed VBS state breaking $\pi/2$ -rotation.

We point out that we find short-range incommensurate spin correlations at $J_{3}=0.45$ in the VBS phase, and the long-range spiral order is established for $J_{3}>0.75$ , as we have discussed in the previous part. The spiral order in the VBS phase is short-ranged as the order parameter is scaled to zero in the thermodynamic limit. As shown in Fig. 7, with increasing system size $L$ , the peak value of spin structure factor at $J_{3}=0.6$ almost stays a constant, indicating the vanishing spin order. In contrast, at $J_{3}=0.8$ the peak value diverges quickly, showing the development of a spin order. In addition, we have also computed the spin correlation function on a $12\times 28$ long strip for further check. In Fig. 8, the spin-spin correlations at $J_{3}=0.35$ (in the gapless QSL phase), $J_{3}=0.45,0.5$ (in the VBS phase), and $J_{3}=0.8$ (in the spiral phase) are shown. In contrast to the power-law decay at $J_{3}=0.35$ (see more results in Ref. Liu et al. (2022b)), the spin correlations at $J_{3}=0.5$ show a clear exponential decay with an oscillation. Note that at $J_{3}=0.45$ and $0.5$ , correlators $(-1)^{r}\langle S_{0}\cdot S_{r}\rangle$ have positive and negative values ¹¹1There is a typo in Fig. 4(a) in Ref. Liu et al. (2022b); the y-axis label should read $|{\langle}{S_{0}}{\cdot}{S_{r}}{\rangle}|$ , rather than ${(-1)^{r}}{\langle}{S_{0}}{\cdot}{S_{r}}{\rangle}$ . This is a signature of the short-range spiral order (although with a very short correlation length) while, in the long-range spiral ordered phase, the spin correlations oscillate but remain finite at long distance.

To visualize the spin pattern, we compute the values of $\langle S_{i_{x},i_{y}}^{x}\rangle$ , $\langle S_{i_{x},i_{y}}^{y}\rangle$ and $\langle S_{i_{x},i_{y}}^{z}\rangle$ on each site at $J_{3}=0.8$ (set $J_{1}=1$ ) and $J_{3}/J_{1}\rightarrow\infty$ (set $J_{1}=0$ ). In Fig. 9(a) and Fig. 10, $\langle S_{i_{x},i_{y}}^{z}\rangle$ on each site are presented. At $J_{3}=0.8$ we see a clear indication of incommensurate long-range spiral order, while for $J_{3}/J_{1}\rightarrow\infty$ it is commensurate with a period of $4$ lattice spacings along both $x$ and $y$ directions and the 3rd-neighbour spin pairs being antiparallel. The spiral orders can be more explicitly visualized by also considering the $x$ and $y$ spin components. We find that the magnitude of the $y$ -components $\langle S^{y}_{i_{x},i_{y}}\rangle$ is extremely small, which is at least three orders weaker than $\langle S^{z}_{i_{x},i_{y}}\rangle$ or $\langle S^{x}_{i_{x},i_{y}}\rangle$ according to our resolution, indicating that the spiral order is nearly coplanar (see Figs. 10(c) and 10(d)). The incommensurate (commensurate) long-range spiral order is also revealed by the magnetic Bragg peaks at two (four) incommensurate (commensurate) wave vectors, as shown in Fig. 5(g) and Fig. 9(b) (Fig. 5(h)).

II.2 $J_{2}\neq 0$

Based on the above results at $J_{2}=0$ , we establish a good understanding for the properties of the VBS phase and spiral order. Now we turn to the $J_{2}\neq 0$ case to explore the whole phase diagram of the spin- $1/2$ $J_{1}$ - $J_{2}$ - $J_{3}$ model. First of all, we focus on two typical cases at $J_{2}=0.3$ and $J_{2}=0.5$ . At $J_{2}=0.3$ , increasing $J_{3}$ from zero, the system sequentially experiences a Néel AFM phase, a gapless QSL phase, and a VBS phase. Note that for $J_{2}=0.5$ there is no Néel AFM phase for $J_{3}\geq 0$ . In both cases, when $J_{3}$ is large, the system always lies in the long-range spiral ordered phase. Here we compute the physical quantities on the $12\times 12$ cluster to check whether there exists other possible phases between the VBS and spiral order phase.

In Fig. 11, we present the variation of the spin and dimer order parameters with increasing $J_{3}$ . We note that at $J_{2}=0.3$ and $J_{3}=0$ (in the Néel AFM phase) the AFM order parameter $\langle S_{z}^{\rm AFM}\rangle$ on the $12\times 12$ system size is still zero, in contrast to the results on the $16\times 16$ system at $J_{2}=0$ . This is because the size $L=12$ is still a bit too small to observe the spontaneous breaking of SU(2) symmetry occurring in the Néel phase. The other spin order parameter $\langle S_{z}^{\rm abs}\rangle$ shows a sharp change around $J_{3}=0.7$ for both $J_{2}=0.3$ and $J_{2}=0.5$ , indicating a first-order transition to the long-range spiral ordered phase. Additionally, the dimer order parameters $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ tend to deviate from each other when $J_{3}$ is large enough, and then have sudden changes at the first-order transitions.

To further investigate the properties of the VBS region, we extend our analysis to larger system sizes, including $L=16$ and $L=20$ . Taking $J_{2}=0.3$ as an example, as shown in Fig.12(a), we observe distinct behaviors at different $J_{3}$ : nonzero $\langle D_{x}\rangle^{2}=\langle D_{x}\rangle^{2}$ in the 2D limit for smaller $J_{3}$ values ( $J_{3}=0.4$ ), and $\langle D_{x}\rangle^{2}\neq\langle D_{y}\rangle^{2}$ for larger $J_{3}$ values ( $J_{3}=0.5$ ), which are consistent with the characteristics of a plaquette VBS and a mixed VBS, respectively. By comparing the two cases $(J_{2},J_{3})=(0.3,0.4)$ and $(J_{2},J_{3})=(0.3,0.5)$ , the transition point between the two VBS states is roughly located at $J_{3}=0.45(5)$ . Similar features in the VBS region are also found for $J_{2}=0.5$ , with a transition point estimated at $J_{3}=0.25(5)$ by comparing the extrapolated thermodynamic limit results of $(J_{2},J_{3})=(0.5,0.2)$ and $(J_{2},J_{3})=(0.5,0.3)$ , shown in Fig.12(b). In fact, the plaquette and mixed VBS state persist across most of the VBS region, and their transition line with respect to $(J_{2},J_{3})$ is represented as the violet line in Fig. 1. We remark that except aforementioned cases, for other values of $(J_{2},J_{3})$ , the transition points (the violet points in Fig. 1) are estimated based on the results of the $16\times 16$ system.

We have also checked the evolution of spin structure factor with respect to $J_{3}$ , as shown in Fig. 13. We find that the maximum of spin structure factor gradually moves away from $(\pi,\pi)$ . A comparison to the results obtained at $J_{2}=0$ reveals the similar features for spin and dimer order parameters as well as for spin structure factor, indicating no other phases between the VBS and spiral order phase. Interestingly, as discussed above, the emergence of short-range spiral spin correlation in some region of the plaquette VBS phase may be considered as a precursor of the transitions to the mixed VBS phase and to the spiral phase (see the transition lines reported in Fig. 1). We point out that, in contrast to the classical case showing the existence of a spiral phase ordered at the wave vector $(\pi,\pm q)$ or $(\pm q,\pi)$ , either our PEPS or previous exact diagonalization results do not support such a long-range ordered spiral phase in the quantum spin- $1/2$ model.

The phase diagram of the $J_{1}$ - $J_{2}$ - $J_{3}$ model is becoming more precise now: we know it must contain the Néel AFM, gapless QSL, VBS, stripe, and long-range spiral order phases. Since the VBS-spiral and VBS-stripe phase transitions are both first order, we can also use the exact ground-state energies on a $4\times 4$ periodic system to evaluate the phase boundaries (see Fig. 18 in Appendix. C). We also use the PEPS results on the larger $12\times 12$ or $16\times 16$ systems with open boundaries to evaluate the VBS-spiral transition points for several $(J_{2},J_{3})$ values, showing agreement with the corresponding estimations from the periodic $4\times 4$ cluster. The estimated phase boundaries can be seen in Fig. 1.

II.3 AFM-spiral transition

It has been already shown in our previous work that for a negative $J_{2}\lesssim-0.25$ , there exists an AFM-VBS transition line Liu et al. (2022b). Here, we find that for a larger negative $J_{2}$ such as $J_{2}=-1.2$ , the VBS phase will disappear and a direct first-order AFM-spiral phase transition occurs.

We first consider the intermediate negative $J_{2}$ values, say $J_{2}=-0.8$ or $J_{2}=-1.0$ . In such cases, increasing $J_{3}$ will lead to a continuous AFM-VBS transition followed by a first-order VBS-spiral transition. As discussed in Sec. II.1.3, the VBS-spiral transition point can be generically determined by comparing the energies of the simulations with different initial states. In Figs. 14(a) and 14(b), we show the energy variation with respect to $J_{3}$ on the $16\times 16$ size when a VBS-spiral transition occurs, at fixed $J_{2}=-0.8$ and $J_{2}=-1.0$ respectively. The AFM-VBS transition point can be evaluated by the finite-size scaling of the corresponding order parameters. In Figs. 14(c) and 14(d), we present the size scaling of the Néel AFM order parameter $\langle M^{2}_{0}\rangle$ and VBS order parameter $\langle D\rangle^{2}=\langle D_{x}\rangle^{2}+\langle D_{y}\rangle^{2}$ for $J_{2}=-1.0$ as an example, where $\langle M^{2}_{0}\rangle=\frac{1}{L^{2}}S({\bf k_{0}})$ with ${\bf k_{0}}=(\pi,\pi)$ . We find that the VBS phase is located in the interval $0.92\lesssim J_{3}\lesssim 0.99$ for $J_{2}=-1.0$ . Upon further increasing the magnitude of the negative $J_{2}$ , the extension of the VBS region with $J_{3}$ shrinks rapidly, leading to an end point of the continuous AFM-VBS transition.

Now we turn to $J_{2}=-1.2$ , at which we find a direct first-order AFM-spiral transition with increasing $J_{3}$ . As shown in Fig. 15, we can see the transition happening at $J_{3}\simeq 1.06$ . At $J_{3}=1.06$ , the optimal AFM state has an energy $E=-0.79285$ , very close to the optimal spiral state energy $E=-0.79342$ (see Fig. 15(a)), but they have clear AFM and spiral spin patterns correspondingly (not shown here). The ground-state spin order parameters also show sharp changes around $J_{3}=1.06$ in Fig. 15(b), providing strong evidence for a first-order phase transition.

II.4 $J_{3}<0$

To obtain more information of the phase diagram, we finally consider the region with a negative $J_{3}$ . For $J_{3}=0$ , i.e. the $J_{1}$ - $J_{2}$ model, a gapless QSL and a VBS phase emerge between the Néel AFM and stripe AFM phases. A negative $J_{3}$ will enhance spin orders including the Néel and stripe AFM orders, which therefore destabilizes the QSL and VBS phases. Here, we try to estimate how large a (negative) $J_{3}$ will be able to suppress the QSL and VBS phases.

We first consider a relatively larger negative $J_{3}$ , $J_{3}=-0.1$ . In this case, we find a direct first-order transition between the Néel AFM and stripe phases. At $J_{2}=0.50$ , we observe a clear Néel AFM pattern and the averaged local moment $\langle S^{z}_{i_{x},i_{y}}\rangle$ is about $0.026$ on the $16\times 16$ cluster. Further increasing $J_{2}$ to $0.58$ , this averaged moment is reduced to $0.016$ , still showing a clear Néel AFM order. Meanwhile, the local dimer order $\langle D\rangle^{2}=\langle D_{x}\rangle^{2}+\langle D_{y}\rangle^{2}$ remains much smaller, which is about $0.0067$ . The Néel-stripe is a typical first-order transition, and the $J_{2}$ transition point is expected to shift as the system size increases, similarly to the VBS-stripe transition in the $J_{1}$ - $J_{2}$ model Liu et al. (2022a). As discussed in Sec. II.1.3, we can evaluate the transition point by initializing the PEPS optimization from either the Néel or the stripe state. Using the crossing of the energy curves shown in Fig. 16(a), the transition points $J_{2}^{c}(L)$ can be obtained as $0.656$ , $0.630$ , and $0.618$ for $8\times 8$ , $12\times 12$ and $16\times 16$ , respectively. Furthermore, we take a linear extrapolation of $J_{2}^{c}(L)$ versus $1/L$ , giving the first-order transition point in the thermodynamic limit $J_{2}^{c}(L\rightarrow\infty)\simeq 0.58$ , as shown in Fig. 16(b). As the bulk energy on the $16\times 16$ size can better approximate the energy in the thermodynamic limit, following Ref. Liu et al. (2022a), we also use the central $4\times 4$ bulk energy of the $16\times 16$ open system to estimate the thermodynamic limit transition point, which gives a consistent value $J_{2}^{c}(L\rightarrow\infty)\simeq 0.582$ (not shown here). We also demonstrate the spin and dimer order parameters on the $16\times 16$ size, as shown in Fig. 16(c), which also confirms the first-order transition nature.

We further consider a smaller $J_{3}$ , say, $J_{3}=-0.05$ . By taking $J_{2}=0.58$ , we find that the spin order $\langle S^{z}_{i_{x},i_{y}}\rangle$ at $(J_{2},J_{3})=(0.58,-0.05)$ is $0.008$ , which is half of the value obtained at $(J_{2},J_{3})=(0.58,-0.1)$ . This result indicates that $(J_{2},J_{3})=(0.58,-0.05)$ could still be in the Néel AFM phase, but rather close to the nonmagnetic regime. These results suggest that both the gapless QSL and VBS phases may disappear at a quite small negative value of $J_{3}$ . In such a narrow region, it is hard to accurately determine the transition lines of the QSL and VBS phases. Therefore, we schematically show in the inset of Fig. 1 the phase boundaries of the QSL and VBS phases for $J_{3}<0$ ending in two multicritical points. Here we assume that the gapless QSL does not touch the stripe phase but remains separated by the VBS phase. We believe that a first-order transition always occurs between the stripe phase and other phases, in that case the AFM and VBS phases. As the two putative multicritical points are very close, another possible scenario is that the Néel AFM, QSL, VBS and stripe phases are directly connected by a quadruple point.

III Conclusion and discussion

In this work, we establish the global phase diagram of the square-lattice $J_{1}$ - $J_{2}$ - $J_{3}$ model using finite PEPS simulations with careful finite-size scaling. We compute spin and dimer order parameters, as well as spin-spin correlation functions. First, we focus on the identification of the VBS phases and find a novel transition from the 4-fold degenerate plaquette VBS phase (beyond the boundary with the previously discovered QSL phase) to a 8-fold degenerate mixed-VBS phase (e.g. in the $J_{1}$ - $J_{3}$ model with increasing $J_{3}$ ). The mixed-VBS phase can be viewed as spontaneously breaking the point group $C_{4}$ symmetry of the plaquette VBS state. In other words, while the horizontal and vertical dimers have the same magnitude in the plaquette VBS phase, they start to become different at the transition to the mixed-VBS. The mixed-VBS phase also shows incommensurate short-range spin correlations while approaching the magnetic spiral phase ordered at wave vectors $(\pm q,\pm q)$ . Our results combined with previous studies suggest that the transition between the plaquette and mixed columnar-plaquette VBS phases is continuous as expected from the Ginzburg-Landau paradigm and studies of quantum dimer models Ralko et al. (2008). In contrast, the transition between the mixed-VBS and spiral phase appears to be first order. In the spiral phase, we establish the existence of long-range spiral spin correlation and explicitly visualize the incommensurate spiral patterns in real space.

Therefore, the overall ground-state phase diagram of the $J_{1}$ - $J_{2}$ - $J_{3}$ model is elucidated in detail, which contains six phases: a Néel AFM phase ordered at $(\pi,\pi)$ , a stripe AFM phase ordered at $(\pi,0)$ or $(0,\pi)$ , two VBS phases, a gapless QSL phase, and a spiral phase ordered at $(\pm q,\pm q)$ including the state ordered at $(\pm\pi/2,\pm\pi/2)$ . We would like to stress that our work significantly broadens the knowledge on the VBS and spiral order phase. It also provides a canonical example for the understanding of quantum effects and competition caused by frustration. Interestingly, we find that the QSL and VBS phases in the $J_{1}$ - $J_{2}$ model, which are easily suppressed by a very small negative $J_{3}=-0.02\sim-0.05$ , are rather close to the multicritical points occurring at a quite small $J_{3}<0$ , as seen in Fig. 1. This closeness may naturally explain the very long correlation lengths found in the previous studies of the pure $J_{1}$ - $J_{2}$ model for the nonmagnetic region Gong et al. (2014); Poilblanc and Mambrini (2017); Liu et al. (2022a).

Furthermore, we find two distinct types of multicritical points in the phase diagram. The first type involves the critical points at which three continuous transition lines intersect: AFM-VBS, AFM-QSL, and QSL-VBS (marked as the two blue dots in Fig. 1). The second type encompasses the critical points that mark the termination of a continuous line culminating in a first-order transition, which occurs with the AFM-VBS critical line reaching either the spiral phase or the stripe phase (denoted as the two red dots in Fig. 1). It is noteworthy that these two types of critical points are conceptually different. Intriguingly, the continuous AFM-VBS transition functioning as a line of deconfined quantum critical points (DQCP) Liu et al. (2022b), can culminate at both types of multicritical points. This observation can offer valuable insights into comprehending the nature of the DQCP Senthil et al. (2004a, b).

Finally, we would like to point out that our study further demonstrates the capability of finite PEPS as a powerful numerical tool to study strongly correlated 2D quantum many-body systems. Since the tensor elements of finite PEPS can be independent, the finite PEPS constitutes a very general ansatz family and can naturally capture nonuniform properties. In particular, the power of finite PEPS has been fully explored when it is used to accurately represent incommensurate short-range or long-range orders. We note that in fermionic correlated systems like the $t-J$ model, short-range incommensurate correlations can also exist Moreo et al. (1990), and hence the finite PEPS method should be able to provide accurate results for such systems.

IV Acknowledgment

We thank Gang Chen, Zheng-Xin Liu, Han-Qiang Wu and Rong Yu for helpful discussions. This work is supported by the CRF C7012-21GF, the ANR/RGC Joint Research Scheme No. A-CUHK402/18 from the Hong Kong’s Research Grants Council and the TNTOP ANR-18-CE30-0026-01 grants awarded from the French Research Council. Wei-Qiang Chen is supported by the National Key R&D Program of China (Grants No. 2022YFA1403700), NSFC (Grants No. 12141402), the Science, Technology and Innovation Commission of Shenzhen Municipality (No. ZDSYS20190902092905285), Guangdong Basic and Applied Basic Research Foundation under Grant No. 2020B1515120100, and Center for Computational Science and Engineering at Southern University of Science and Technology. S.S.G. was supported by the NSFC (No. 12274014), the Special Project in Key Areas for Universities in Guangdong Province (No. 2023ZDZX3054), and the Dongguan Key Laboratory of Artificial Intelligence Design for Advanced Materials (DKL-AIDAM). W.Y.L. was supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator.

Appendix A convergence with $D$

Here we check the $D$ -convergence on a $16\times 16$ open system at $(J_{2},J_{3})=(-1,0.8)$ . In Table.1, we present the obtained $D=4-10$ results for energy, dimerization and magnetization. It clearly indicates $D=8$ is sufficient to get converged results.

Table 1: The convergence of physical quantities with respect to the PEPS bond dimension

D

on a

16\times 16

system at

(J_{2},J_{3})=(-1,0.8)

, including ground state energy per site, dimerization

\langle D\rangle^{2}=\langle D_{x}\rangle^{2}+\langle D_{y}\rangle^{2}

, and AFM magnetization

\langle M^{2}_{z}\rangle=\langle(S^{\rm AFM}_{z})^{2}\rangle

(see main text for detailed definitions).

$D$	$E$	$\langle D\rangle^{2}$	$\langle M^{2}_{z}\rangle$
4	-0.769042(6)	0.00656(9)	0.0290(1)
6	-0.772245(4)	0.00686(8)	0.0261(2)
8	-0.773084(3)	0.00701(6)	0.0217(1)
10	-0.773092(5)	0.00701(8)	0.0218(2)

Appendix B Finite size scaling of VBS order parameters

Figure 17 presents the $1/L$ dependence of the dimer order parameters $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ for the $J_{1}-J_{3}$ model (i.e. $J_{2}$ =0), with different $J_{3}$ values inside the VBS phase. We can see $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ are almost identical on each size for smaller $J_{3}$ , and their extrapolated values for 2D limit are 0.0066(5) and 0.0059(4) for $J_{3}=0.45$ , 0.0122(8) and 0.0121(7) for $J_{3}=0.5$ , correspondingly. This indicates $J_{3}=0.45$ and $J_{3}=0.5$ have a plaquette VBS. While for $J_{3}=0.6$ and 0.7, $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ get gradually different with system size $L$ increasing. Furthermore, The extrapolated values of $\langle D_{x}\rangle^{2}$ and $\langle D_{y}\rangle^{2}$ for 2D limit both are nonzero, 0.0178(16) and 0.0270(7) for $J_{3}=0.6$ , and 0.0276(21) and 0.0378(9) for $J_{3}=0.7$ , consistent with a mixed VBS phase.

Appendix C $4\times 4$ results from exact diagonalization

As the VBS-spiral and VBS-stripe phase transitions are first-order, we use the exact ground state energies on a $4\times 4$ periodic system to estimate the phase boundaries, as shown in Fig. 18. For the VBS-stripe transition, given $J_{3}=0$ , we have found that the transition point is located at $J_{2}=0.61$ Liu et al. (2022a). Here the second derivative of the $4\times 4$ energy with respect to $J_{2}$ gives an estimate of the transition point around $J_{2}=0.626$ , in good agreement. For the VBS-spiral transition at fixed $J_{2}=0$ , 0.3 and 0.5, the transition points obtained by PEPS calculations on open $16\times 16$ or $12\times 12$ systems, are $J_{3}\simeq 0.75$ , 0.70 and 0.68, respectively, also consistent with the corresponding estimations from the periodic $4\times 4$ cluster, namely $J_{3}\simeq 0.757$ , 0.691 and 0.674.

Appendix D spin structure factor in mixed-VBS phase

We have shown the contour plots of the spin structure factor $S(\bf k)$ of the $J_{1}$ - $J_{3}$ model in Fig. 5. Here we have a much closer look at $S(\bf k)$ in the mixed-VBS phase searching for (weak) signatures of the spontaneously breaking of the $\pi/2$ rotation symmetry. In Fig. 19, we present $S(\bf k)$ along $k_{x}=k_{y}$ and $k_{x}=-k_{y}$ at $J_{3}=0.5$ (inside the plaquette VBS phase) and $J_{3}=0.7$ (inside the mixed columnar-plaquette VBS phase). At $J_{3}=0.5$ , the two curves are almost identical, consistent with a $\pi/2$ rotation symmetry for the plaquette VBS. In contrast, at $J_{3}=0.7$ , they show visible differences around the maxima, indicating the breaking of $\pi/2$ rotation symmetry in the mixed VBS phase. Note that a clear incommensurability and a larger spin correlation length (estimated from the inverse of the width of the peaks) are seen in this case.

References

Anderson (1987) P. W. Anderson, “The resonating valence bond state in ${{\rm La_{2}CuO_{4}}}$ and superconductivity,” Science 235, 1196–1198 (1987).
Gelfand et al. (1989) Martin P. Gelfand, Rajiv R. P. Singh, and David A. Huse, “Zero-temperature ordering in two-dimensional frustrated quantum Heisenberg antiferromagnets,” Phys. Rev. B 40, 10801–10809 (1989).
Moreo et al. (1990) Adriana Moreo, Elbio Dagotto, Thierry Jolicoeur, and José Riera, “Incommensurate correlations in the $t$ - ${J}$ and frustrated spin-1/2 Heisenberg models,” Phys. Rev. B 42, 6283–6293 (1990).
Chandra and Douçot (1988) P. Chandra and B. Douçot, “Possible spin-liquid state at large ${S}$ for the frustrated square Heisenberg lattice,” Phys. Rev. B 38, 9335–9338 (1988).
Figueirido et al. (1990) F. Figueirido, A. Karlhede, S. Kivelson, S. Sondhi, M. Rocek, and D. S. Rokhsar, “Exact diagonalization of finite frustrated spin-1/2 Heisenberg models,” Phys. Rev. B 41, 4619–4632 (1990).
Chubukov (1991) Andrey Chubukov, “First-order transition in frustrated quantum antiferromagnets,” Phys. Rev. B 44, 392–394 (1991).
Read and Sachdev (1991) N. Read and Subir Sachdev, “Large-N expansion for frustrated quantum antiferromagnets,” Phys. Rev. Lett. 66, 1773–1776 (1991).
Ioffe and Larkin (1988) L. B. Ioffe and A. I. Larkin, “Effective action of a two-dimensional antiferromagnet,” Int. J. Mod. Phys. B 2, 203–219 (1988).
Einarsson and Johannesson (1991) Torbjörn Einarsson and Henrik Johannesson, “Effective-action approach to the frustrated Heisenberg antiferromagnet in two dimensions,” Phys. Rev. B 43, 5867–5882 (1991).
Ferrer (1993) Jaime Ferrer, “Spin-liquid phase for the frustrated quantum Heisenberg antiferromagnet on a square lattice,” Phys. Rev. B 47, 8769–8782 (1993).
Capriotti and Sachdev (2004) Luca Capriotti and Subir Sachdev, “Low-temperature broken-symmetry phases of spiral antiferromagnets,” Phys. Rev. Lett. 93, 257206 (2004).
Leung and Lam (1996) P. W. Leung and Ngar-wing Lam, “Numerical evidence for the spin-Peierls state in the frustrated quantum antiferromagnet,” Phys. Rev. B 53, 2213–2216 (1996).
Capriotti et al. (2004) Luca Capriotti, Douglas J. Scalapino, and Steven R. White, “Spin-liquid versus dimerized ground states in a frustrated Heisenberg antiferromagnet,” Phys. Rev. Lett. 93, 177004 (2004).
Mambrini et al. (2006) Matthieu Mambrini, Andreas Läuchli, Didier Poilblanc, and Frédéric Mila, “Plaquette valence-bond crystal in the frustrated Heisenberg quantum antiferromagnet on the square lattice,” Phys. Rev. B 74, 144422 (2006).
Ralko et al. (2009) A. Ralko, M. Mambrini, and D. Poilblanc, “Generalized quantum dimer model applied to the frustrated Heisenberg model on the square lattice: Emergence of a mixed columnar-plaquette phase,” Phys. Rev. B 80, 184427 (2009).
du Croo de Jongh et al. (2000) M. S. L. du Croo de Jongh, J. M. J. van Leeuwen, and W. van Saarloos, “Incorporation of density-matrix wave functions in monte carlo simulations: Application to the frustrated Heisenberg model,” Phys. Rev. B 62, 14844–14854 (2000).
Sindzingre et al. (2010) Philippe Sindzingre, Nic Shannon, and Tsutomu Momoi, “Phase diagram of the spin-1/2 ${J}_{1}$ - ${J}_{2}$ - ${J}_{3}$ Heisenberg model on the square lattice,” Journal of Physics: Conference Series 200, 022058 (2010).
Dagotto and Moreo (1989) Elbio Dagotto and Adriana Moreo, “Phase diagram of the frustrated spin-1/2 Heisenberg antiferromagnet in 2 dimensions,” Phys. Rev. Lett. 63, 2148–2151 (1989).
Sachdev and Bhatt (1990) Subir Sachdev and R. N. Bhatt, “Bond-operator representation of quantum spins: Mean-field theory of frustrated quantum Heisenberg antiferromagnets,” Phys. Rev. B 41, 9323–9329 (1990).
Poilblanc et al. (1991) Didier Poilblanc, Eduardo Gagliano, Silvia Bacci, and Elbio Dagotto, “Static and dynamical correlations in a spin-1/2 frustrated antiferromagnet,” Phys. Rev. B 43, 10970–10983 (1991).
Chubukov and Jolicoeur (1991) Andrey V. Chubukov and Th. Jolicoeur, “Dimer stability region in a frustrated quantum Heisenberg antiferromagnet,” Phys. Rev. B 44, 12050–12053 (1991).
Schulz and Ziman (1992) H. J Schulz and T. A. L Ziman, “Finite-size scaling for the two-dimensional frustrated quantum Heisenberg antiferromagnet,” Europhysics Letters 18, 355–360 (1992).
Ivanov and Ivanov (1992) N. B. Ivanov and P. Ch. Ivanov, “Frustrated two-dimensional quantum Heisenberg antiferromagnet at low temperatures,” Phys. Rev. B 46, 8206–8213 (1992).
Einarsson and Schulz (1995) T. Einarsson and H. J. Schulz, “Direct calculation of the spin stiffness in the ${J}_{1}\text{$-$}{J}_{2}$ Heisenberg antiferromagnet,” Phys. Rev. B 51, 6151–6154 (1995).
H. J. Schulz and Poilblanc (1996) T. A. L. Ziman H. J. Schulz and D. Poilblanc, “Magnetic order and disorder in the frustrated quantum Heisenberg antiferromagnet in two dimensions,” J. Phys. I 6, 675–703 (1996).
Zhitomirsky and Ueda (1996) M. E. Zhitomirsky and Kazuo Ueda, “Valence-bond crystal phase of a frustrated spin-1/2 square-lattice antiferromagnet,” Phys. Rev. B 54, 9007–9010 (1996).
Singh et al. (1999) Rajiv R. P. Singh, Zheng Weihong, C. J. Hamer, and J. Oitmaa, “Dimer order with striped correlations in the ${J}_{1}$ - ${J}_{2}$ Heisenberg model,” Phys. Rev. B 60, 7278–7283 (1999).
Capriotti and Sorella (2000) Luca Capriotti and Sandro Sorella, “Spontaneous plaquette dimerization in the ${J}_{1}$ - ${J}_{2}$ Heisenberg model,” Phys. Rev. Lett. 84, 3173–3176 (2000).
Capriotti et al. (2001) Luca Capriotti, Federico Becca, Alberto Parola, and Sandro Sorella, “Resonating valence bond wave functions for strongly frustrated spin systems,” Phys. Rev. Lett. 87, 097201 (2001).
Zhang et al. (2003) Guang-Ming Zhang, Hui Hu, and Lu Yu, “Valence-bond spin-liquid state in two-dimensional frustrated spin- $1/2$ Heisenberg antiferromagnets,” Phys. Rev. Lett. 91, 067201 (2003).
Takano et al. (2003) Ken’ichi Takano, Yoshiya Kito, Yoshiaki Ōno, and Kazuhiro Sano, “Nonlinear $\sigma$ model method for the ${J}_{1}$ - ${J}_{2}$ Heisenberg model: Disordered ground state with plaquette symmetry,” Phys. Rev. Lett. 91, 197202 (2003).
Sirker et al. (2006) J. Sirker, Zheng Weihong, O. P. Sushkov, and J. Oitmaa, “ ${J}_{1}$ - ${J}_{2}$ model: First-order phase transition versus deconfinement of spinons,” Phys. Rev. B 73, 184420 (2006).
Schmalfuß et al. (2006) D. Schmalfuß, R. Darradi, J. Richter, J. Schulenburg, and D. Ihle, “Quantum ${J}_{1}$ - ${J}_{2}$ antiferromagnet on a stacked square lattice: Influence of the interlayer coupling on the ground-state magnetic ordering,” Phys. Rev. Lett. 97, 157201 (2006).
Darradi et al. (2008) R. Darradi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Krüger, and J. Richter, “Ground state phases of the spin-1/2 ${J}_{1}$ - ${J}_{2}$ Heisenberg antiferromagnet on the square lattice: A high-order coupled cluster treatment,” Phys. Rev. B 78, 214415 (2008).
Arlego and Brenig (2008) Marcelo Arlego and Wolfram Brenig, “Plaquette order in the ${J}_{1}$ - ${J}_{2}$ - ${J}_{3}$ model: Series expansion analysis,” Phys. Rev. B 78, 224415 (2008).
Isaev et al. (2009) L. Isaev, G. Ortiz, and J. Dukelsky, “Hierarchical mean-field approach to the ${J}_{1}$ - ${J}_{2}$ Heisenberg model on a square lattice,” Phys. Rev. B 79, 024409 (2009).
Murg et al. (2009) V. Murg, F. Verstraete, and J. I. Cirac, “Exploring frustrated spin systems using projected entangled pair states,” Phys. Rev. B 79, 195119 (2009).
Beach (2009) K. S. D. Beach, “Master equation approach to computing rvb bond amplitudes,” Phys. Rev. B 79, 224431 (2009).
Richter and Schulenburg (2010) J. Richter and J. Schulenburg, “The spin-1/2 ${J}_{1}\text{$-$}{J}_{2}$ Heisenberg antiferromagnet on the square lattice: Exact diagonalization for ${{\rm N=40}}$ spins,” The European Physical Journal B 73, 117 (2010).
Yu and Kao (2012) Ji-Feng Yu and Ying-Jer Kao, “Spin- $\frac{1}{2}$ ${J}_{1}$ - ${J}_{2}$ Heisenberg antiferromagnet on a square lattice: A plaquette renormalized tensor network study,” Phys. Rev. B 85, 094407 (2012).
Jiang et al. (2012) Hong-Chen Jiang, Hong Yao, and Leon Balents, “Spin liquid ground state of the spin- $\frac{1}{2}$ square ${J}_{1}$ - ${J}_{2}$ Heisenberg model,” Phys. Rev. B 86, 024424 (2012).
Mezzacapo (2012) Fabio Mezzacapo, “Ground-state phase diagram of the quantum ${J}_{1}$ - ${J}_{2}$ model on the square lattice,” Phys. Rev. B 86, 045115 (2012).
Wang et al. (2013) Ling Wang, Didier Poilblanc, Zheng-Cheng Gu, Xiao-Gang Wen, and Frank Verstraete, “Constructing a gapless spin-liquid state for the spin- $1/2$ ${J}_{1}$ - ${J}_{2}$ Heisenberg model on a square lattice,” Phys. Rev. Lett. 111, 037202 (2013).
Hu et al. (2013) Wen-Jun Hu, Federico Becca, Alberto Parola, and Sandro Sorella, “Direct evidence for a gapless ${Z}_{2}$ spin liquid by frustrating néel antiferromagnetism,” Phys. Rev. B 88, 060402 (2013).
Doretto (2014) R. L. Doretto, “Plaquette valence-bond solid in the square-lattice ${J}_{1}$ - ${J}_{2}$ antiferromagnet Heisenberg model: A bond operator approach,” Phys. Rev. B 89, 104415 (2014).
Qi and Gu (2014) Yang Qi and Zheng-Cheng Gu, “Continuous phase transition from Néel state to ${Z}_{2}$ spin-liquid state on a square lattice,” Phys. Rev. B 89, 235122 (2014).
Gong et al. (2014) Shou-Shu Gong, Wei Zhu, D. N. Sheng, Olexei I. Motrunich, and Matthew P. A. Fisher, “Plaquette ordered phase and quantum phase diagram in the spin- $\frac{1}{2}$ ${J}_{1}$ - ${J}_{2}$ square Heisenberg model,” Phys. Rev. Lett. 113, 027201 (2014).
Chou and Chen (2014) Chung-Pin Chou and Hong-Yi Chen, “Simulating a two-dimensional frustrated spin system with fermionic resonating-valence-bond states,” Phys. Rev. B 90, 041106 (2014).
Morita et al. (2015) Satoshi Morita, Ryui Kaneko, and Masatoshi Imada, “Quantum spin liquid in spin 1/2 ${J}_{1}$ - ${J}_{2}$ Heisenberg model on square lattice: Many-variable variational monte carlo study combined with quantum-number projections,” Journal of the Physical Society of Japan 84, 024720 (2015).
Richter et al. (2015) Johannes Richter, Ronald Zinke, and Damian J. J. Farnell, “The spin-1/2 square-lattice ${J}_{1}$ - ${J}_{2}$ model: the spin-gap issue,” The European Physical Journal B 88, 2 (2015).
Wang et al. (2016) Ling Wang, Zheng-Cheng Gu, Frank Verstraete, and Xiao-Gang Wen, “Tensor-product state approach to spin- $\frac{1}{2}$ square ${J}_{1}$ - ${J}_{2}$ antiferromagnetic Heisenberg model: Evidence for deconfined quantum criticality,” Phys. Rev. B 94, 075143 (2016).
Poilblanc and Mambrini (2017) Didier Poilblanc and Matthieu Mambrini, “Quantum critical phase with infinite projected entangled paired states,” Phys. Rev. B 96, 014414 (2017).
Wang and Sandvik (2018) Ling Wang and Anders W. Sandvik, “Critical level crossings and gapless spin liquid in the square-lattice spin- $1/2$ ${J}_{1}$ - ${J}_{2}$ Heisenberg antiferromagnet,” Phys. Rev. Lett. 121, 107202 (2018).
Haghshenas and Sheng (2018) R. Haghshenas and D. N. Sheng, “ $u(1)$ -symmetric infinite projected entangled-pair states study of the spin-1/2 square ${J}_{1}$ - ${J}_{2}$ Heisenberg model,” Phys. Rev. B 97, 174408 (2018).
Liu et al. (2018) Wen-Yuan Liu, Shaojun Dong, Chao Wang, Yongjian Han, Hong An, Guang-Can Guo, and Lixin He, “Gapless spin liquid ground state of the spin-1/2 ${J}_{1}$ - ${J}_{2}$ Heisenberg model on square lattices,” Phys. Rev. B 98, 241109 (2018).
Poilblanc et al. (2019) Didier Poilblanc, Matthieu Mambrini, and Sylvain Capponi, “Critical colored-RVB states in the frustrated quantum Heisenberg model on the square lattice,” SciPost Phys. 7, 41 (2019).
Hasik et al. (2021) Juraj Hasik, Didier Poilblanc, and Federico Becca, “Investigation of the néel phase of the frustrated Heisenberg antiferromagnet by differentiable symmetric tensor networks,” SciPost Phys. 10, 12 (2021).
Ferrari and Becca (2020) Francesco Ferrari and Federico Becca, “Gapless spin liquid and valence-bond solid in the ${J}_{1}$ - ${J}_{2}$ Heisenberg model on the square lattice: Insights from singlet and triplet excitations,” Phys. Rev. B 102, 014417 (2020).
Liu et al. (2022a) Wen-Yuan Liu, Shou-Shu Gong, Yu-Bin Li, Didier Poilblanc, Wei-Qiang Chen, and Zheng-Cheng Gu, “Gapless quantum spin liquid and global phase diagram of the spin-1/2 ${J}_{1}$ - ${J}_{2}$ square antiferromagnetic Heisenberg model,” Science Bulletin 67, 1034–1041 (2022a).
Nomura and Imada (2021) Yusuke Nomura and Masatoshi Imada, “Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy,” Phys. Rev. X 11, 031034 (2021).
Senthil et al. (2004a) T. Senthil, Ashvin Vishwanath, Leon Balents, Subir Sachdev, and Matthew P. A. Fisher, “Deconfined quantum critical points,” Science 303, 1490–1494 (2004a).
Liu et al. (2022b) Wen-Yuan Liu, Juraj Hasik, Shou-Shu Gong, Didier Poilblanc, Wei-Qiang Chen, and Zheng-Cheng Gu, “Emergence of gapless quantum spin liquid from deconfined quantum critical point,” Phys. Rev. X 12, 031039 (2022b).
Wu et al. (2022) Muwei Wu, Shou-Shu Gong, Dao-Xin Yao, and Han-Qing Wu, “Phase diagram and magnetic excitations of ${J}_{1}$ - ${J}_{3}$ Heisenberg model on the square lattice,” Phys. Rev. B 106, 125129 (2022).
Corboz et al. (2011) Philippe Corboz, Steven R. White, Guifré Vidal, and Matthias Troyer, “Stripes in the two-dimensional $t$ - ${J}$ model with infinite projected entangled-pair states,” Phys. Rev. B 84, 041108 (2011).
Corboz et al. (2014) Philippe Corboz, T. M. Rice, and Matthias Troyer, “Competing states in the $t$ - ${J}$ model: Uniform d-wave state versus stripe state,” Phys. Rev. Lett. 113, 046402 (2014).
Zheng et al. (2017) Bo-Xiao Zheng, Chia-Min Chung, Philippe Corboz, Georg Ehlers, Ming-Pu Qin, Reinhard M Noack, Hao Shi, Steven R White, Shiwei Zhang, and Garnet Kin-Lic Chan, “Stripe order in the underdoped region of the two-dimensional Hubbard model,” Science 358, 1155–1160 (2017).
Liu et al. (2017) Wen-Yuan Liu, Shao-Jun Dong, Yong-Jian Han, Guang-Can Guo, and Lixin He, “Gradient optimization of finite projected entangled pair states,” Phys. Rev. B 95, 195154 (2017).
Liu et al. (2021) Wen-Yuan Liu, Yi-Zhen Huang, Shou-Shu Gong, and Zheng-Cheng Gu, “Accurate simulation for finite projected entangled pair states in two dimensions,” Phys. Rev. B 103, 235155 (2021).
Liu et al. (2022c) Wen-Yuan Liu, Shou-Shu Gong, Wei-Qiang Chen, and Zheng-Cheng Gu, “Emergent symmetry in frustrated magnets: From deconfined quantum critical point to gapless quantum spin liquid,” arXiv:2212.00707 (2022c).
Liu et al. (2023) Wen-Yuan Liu, Xiao-Tian Zhang, Zhe Wang, Shou-Shu Gong, Wei-Qiang Chen, and Zheng-Cheng Gu, “Deconfined quantum criticality with emergent symmetry in the extended shastry-sutherland model,” arXiv:2309.10955 (2023).
Ralko et al. (2008) A. Ralko, D. Poilblanc, and R. Moessner, “Generic mixed columnar-plaquette phases in Rokhsar-Kivelson models,” Phys. Rev. Lett. 100, 037201 (2008).
Jiang et al. (2008) H. C. Jiang, Z. Y. Weng, and T. Xiang, “Accurate determination of tensor network state of quantum lattice models in two dimensions,” Phys. Rev. Lett. 101, 090603 (2008).
Locher (1990) Peter Locher, “Linear spin waves in a frustrated Heisenberg model,” Phys. Rev. B 41, 2537–2539 (1990).
Note (1) There is a typo in Fig. 4(a) in Ref. Liu et al. (2022b); the y-axis label should read $|{\langle}{S_{0}}{\cdot}{S_{r}}{\rangle}|$ , rather than ${(-1)^{r}}{\langle}{S_{0}}{\cdot}{S_{r}}{\rangle}$ .
Senthil et al. (2004b) T. Senthil, Leon Balents, Subir Sachdev, Ashvin Vishwanath, and Matthew P. A. Fisher, “Quantum criticality beyond the Landau-Ginzburg-Wilson paradigm,” Phys. Rev. B 70, 144407 (2004b).

Tensor network study of the spin-1/2 square-lattice J1J_{1}-J2J_{2}-J3J_{3} model: incommensurate spiral order, mixed valence-bond solids, and multicritical points