Is the Valence Bond Solid state in $J_{1}$-$J_{2}$ Square Lattice Heisenberg Model Plaquette or Columnar?

Abstract

We utilize Density Matrix Renormalization Group (DMRG) and Fully Augmented Matrix Product States (FAMPS) methods to investigate the Valence Bond Solid (VBS) phase in the $J_{1}$-$J_{2}$ square lattice Heisenberg model. To differentiate between the Columnar Valence Bond Solid (CVBS) and Plaquette Valence Bond Solid (PVBS) phases, we introduce an anisotropy $\Delta_{y}$ in the nearest neighboring coupling in the $y$-direction, aiming at detecting the possible spontaneous rotational symmetry breaking in the VBS phase. In the calculations, we push the bond dimension to as large as $D=25000$ in FAMPS, simulating systems at a maximum size of $14\times 14$. With a careful extrapolation of the truncation errors and appropriate finite-size scaling, followed by finite $\Delta_{y}$ scaling analysis of the VBS dimer order parameters, we identify the VBS phase as a PVBS type, meaning there is no spontaneous rotational symmetry breaking in the VBS phase. This study not only resolves the long-standing issue of the characterization of the VBS order in the $J_{1}$-$J_{2}$ square lattice Heisenberg model but also highlights the capabilities of FAMPS in the study of two-dimensional quantum many-body systems.

I Introduction

The $J_{1}$-$J_{2}$ spin-$1/2$ square lattice Heisenberg model stands as one of the most extensively researched paradigms for studying the frustration effect in quantum many-body systems. The competing antiferromagnetic interactions between the nearest and next-nearest neighbors give rise to a rich variety of phases, which makes this model a well-known playground for searching exotic quantum states, such as Quantum Spin Liquid (QSL) and Valence Bond Solid (VBS) Dagotto and Moreo (1989); Gong et al. (2014); Capriotti et al. (2001); Read and Sachdev (1991). Understanding these exotic states may be crucial for elucidating the underlying physics of high-temperature superconductors and other strongly correlated materials Anderson (1987); Wen et al. (1989); Read and Sachdev (1991); Lee et al. (2006).

Over the past few decades, extensive research has been conducted to investigate the phase diagram of the $J_{1}$-$J_{2}$ square lattice Heisenberg model. A general agreement has been established that when $J_{2}\rightarrow 0$, the model’s ground state manifests as a Néel antiferromagnetic order Sandvik (1997), which extends to a finite region of $J_{2}/J_{1}\approx 0.5$. For very large but finite values of $J_{2}/J_{1}$, the model’s ground state has antiferromagnetic (AFM) stripe order Dagotto and Moreo (1989); Morita et al. (2015). Intriguingly, the intermediate regime of $0.5\lesssim J_{2}/J_{1}\lesssim 0.6$ is termed the non-magnetic phase, which has been a focal point of ongoing research. Many methods have been employed to study this model, including the exact diagonalization Figueirido et al. (1990); Dagotto and Moreo (1989); Poilblanc et al. (1991); H.J. Schulz et al. (1996); Mambrini et al. (2006a), series expansion Oitmaa and Weihong (1996); Singh et al. (1999); Sirker et al. (2006), quantum Monte Carlo Morita et al. (2015); Hu et al. (2013); Ferrari and Becca (2020), and tensor network methods Jiang et al. (2012); Gong et al. (2014); Wang and Sandvik (2018); Haghshenas and Sheng (2018); Hasik et al. (2021); Liu et al. (2022); Nomura and Imada (2021); Poilblanc and Mambrini (2017).

Despite the valuable insights gained from these results, which have significantly enhanced our understanding of the phase diagram of the $J_{1}$-$J_{2}$ square lattice Heisenberg model, the nature of the non-magnetic regime continues to be a subject of active discussion. A variety of competing states have been suggested for this regime with different methods. These potential states encompass the plaquette valence bond solid (PVBS) Capriotti and Sorella (2000); Takano et al. (2003); Mambrini et al. (2006a, b); Darradi et al. (2008); Yu and Kao (2012); Doretto (2014); Gong et al. (2014); Ferrari and Becca (2020); Nomura and Imada (2021); Wang and Sandvik (2018), columnar valence bond solid (CVBS) Dagotto and Moreo (1989); Poilblanc et al. (1991); Schulz and Ziman (1992); Haghshenas and Sheng (2018); Singh et al. (1999); Read and Sachdev (1991), and quantum spin liquid states Capriotti et al. (2001); Zhang et al. (2003); Wang et al. (2013); Hu et al. (2013); Morita et al. (2015); Poilblanc and Mambrini (2017); Wang and Sandvik (2018); Liu et al. (2018); Ferrari and Becca (2020); Nomura and Imada (2021). The source of these controversies can be attributed to the difficulty of accurately simulating large-scale quantum many-body systems.

Recently, a novel approach, termed Fully Augmented Matrix Product States (FAMPS), has been developed to address and mitigate the limitations inherent in the Density Matrix Renormalization Group (DMRG) method when applied to the study of two-dimensional systems Qian and Qin (2023). This method has been subsequently applied to investigate the phase diagram of $J_{1}$-$J_{2}$ square lattice Heisenberg model Qian and Qin (2024). In Qian and Qin (2024), it is found that a VBS phase directly connects Néel and stripe AFM phases, indicating the absence of a spin liquid phase in the $J_{1}$-$J_{2}$ square lattice Heisenberg model (see Fig. 1 (b)). However, the exact nature of the VBS phase, whether it is a PVBS or a CVBS, remains undetermined. To address this long-standing issue, we employ the DMRG and FAMPS methods, to investigate the nature of the VBS phase at $J_{2}=0.57$, which is deep in the VBS phase.

By applying an anisotropy $\Delta_{y}$ in the nearest neighboring coupling in the $y$-direction in the $J_{1}$-$J_{2}$ square lattice Heisenberg model, we can detect the possible spontaneous rotational symmetry breaking of the VBS phase, which can help us to distinguish CVBS from PVBS. We push the bond dimension to $D=25000$ in our FAMPS calculations and simulate systems with a maximum size of $14\times 14$. Through meticulous extrapolation with truncation errors and reliable finite-size and finite $\Delta_{y}$ scaling analysis, we establish that the VBS phase in the $J_{1}$-$J_{2}$ square lattice Heisenberg model is a PVBS type.

The rest of this paper is structured as follows: in Sec. II, we introduce the model in which we add an anisotropy $\Delta_{y}$ in the nearest neighboring coupling in the $y$-direction to the $J_{1}$-$J_{2}$ square lattice Heisenberg model. We also outline the methods used to study this model and demonstrate the advantage of the FAMPS method in comparison to the pure DMRG method. In Sec. III, we present the results of the VBS dimer order parameter and discuss the approach to determine the nature of the VBS phase. We conclude our work in Sec. IV.

II Model and Methods

II.1 Model

The Hamiltonian of $J_{1}$-$J_{2}$ square lattice Heisenberg model is given by

$$H=J_{1}\sum_{\langle i,j\rangle}\hat{S}_{i}\cdot\hat{S}_{j}+J_{2}\sum_{\langle\langle i,j\rangle\rangle}\hat{S}_{i}\cdot\hat{S}_{j}$$ (1)

with $\langle i,j\rangle$ denoting the nearest neighboring sites and $\langle\langle i,j\rangle\rangle$ denoting the next-nearest neighboring sites. $\hat{S_{i}}$ is the spin-1/2 operator at site $i$ and $J_{1}$ and $J_{2}$ are the nearest and next-nearest neighboring couplings, respectively.

Our study focuses on a square lattice with $L_{x}=L_{y}$ with open boundary conditions in both directions, which preserves the rotational symmetry. The ground state of the VBS phase could potentially manifest as a twofold degenerate CVBS or a genuine PVBS. It is known that there is no spontaneous symmetry breaking for systems with finite size, so it is hard to distinguish between these two states directly in finite systems. To address this issue, we introduce an anisotropy $\Delta_{y}$ in the $y$-direction of the nearest neighboring coupling to detect the possible spontaneous rotational symmetry breaking of the VBS state in the thermodynamic limit. By first extrapolating the system sizes to the thermodynamic limit, followed by the extrapolation to the zero $\Delta_{y}$ limit, we can determine the true nature of the VBS phase. The new Hamiltonian is given by

$$\begin{split}H=\sum_{r}&J_{1}\hat{S}_{r}\cdot\hat{S}_{r+\hat{x}}+(J_{1}+\Delta_{y})\hat{S}_{r}\cdot\hat{S}_{r+\hat{y}}\\ &+J_{2}(\hat{S}_{r}\cdot\hat{S}_{r+\hat{x}+\hat{y}}+\hat{S}_{r}\cdot\hat{S}_{r+\hat{x}-\hat{y}})\end{split}$$ (2)

When $\Delta_{y}$ is set to $0$, this model reduces to the $J_{1}$-$J_{2}$ square lattice Heisenberg model. In this work, we set $J_{1}=1.0$ as the energy unit and $J_{2}=0.57$, deep in the VBS phase as shown in Fig. 1(b).

Figure 1: (a) Illustration of an $8\times 8$ $J_{1}$-$J_{2}$ Heisenberg model on a square lattice, where black lines indicate nearest neighboring couplings $J_{1}$ and dashed blue lines represent next-nearest neighboring couplings $J_{2}$. Open boundary conditions are applied in both $x$ and $y$ directions. The central zone marked by a red dashed area highlights the lattice region utilized for calculating the VBS dimer order parameter $D_{\alpha}$ defined in Eq. (5) (for each system, we use only the results in the central half for the calculation of $D_{\alpha}$). (b) Phase diagram of the $J_{1}$-$J_{2}$ square lattice Heisenberg model, illustrating a VBS phase sandwiched by the Néel and stripe AFM phases within the coupling ratio range $0.535\lesssim J_{2}/J_{1}\lesssim 0.610$ Qian and Qin (2024). The orange arrow at $J_{2}=0.57$ denotes the specific coupling strength at which the VBS dimer order parameter is calculated in this work.

II.2 Method

DMRG is a powerful numerical method for studying one-dimensional and quasi-one-dimensional quantum systems White (1992, 1993); Schollwöck (2005). It is based on the Matrix Product States (MPS) ansatz which is a variational wavefunction ansatz for one-dimensional quantum systems Östlund and Rommer (1995); Verstraete and Cirac (2006); F. Verstraete and Cirac (2008). The MPS ansatz is defined as

$$|\mathrm{MPS}\rangle=\sum_{\{s_{i}\}}\text{Tr}(A^{s_{1}}A^{s_{2}}\cdots A^{s_{N}})|s_{1}s_{2}\cdots s_{N}\rangle$$ (3)

where $A^{s_{i}}$ is the rank-3 local tensor at site $i$ with one physical index $s_{i}$ (with dimension $d$) and two auxiliary indices (with bond dimension $D$). $D$ is the key parameter in DMRG calculations which determines the accuracy of the calculation. Even though DMRG is designed for one-dimensional systems, it can be extended to two-dimensional systems. Nonetheless, the bond dimension required for the accurate simulation of two-dimensional systems increases exponentially with the system size, due to the area-law of entanglement entropy Eisert et al. (2010). Hence for DMRG calculations, achieving a sufficiently large bond dimension to study a two-dimensional system with high precision becomes computationally challenging.

Recently, FAMPS method has been introduced to address the limitations of DMRG in two-dimensional systems. This approach enhances the entanglement in MPS and improves computational accuracy by incorporating an additional layer of unitary tensors, known as disentanglers Vidal (2008), to the physical indices of the MPS. The FAMPS ansatz is defined as

$$|\mathrm{FAMPS}\rangle=D(u)|\mathrm{MPS}\rangle$$ (4)

where $D(u)$ denotes the additional disentangler layer. FAMPS method has been demonstrated to attain more precise results than pure DMRG method and can support area-law-like entanglement for 2D systems while maintaining the low cost of DMRG $\mathrm{O}(D^{3})$ with a small overhead $\mathrm{O}(d^{4})$ Qian and Qin (2023); Felser et al. (2021); Qian and Qin (2022). In this work, we employ the FAMPS method to study the $J_{1}$-$J_{2}$ square lattice Heisenberg model at $J_{2}=0.57$ with the largest system size $14\times 14$ and largest bond dimension $D=25000$. The other system sizes are studied with pure DMRG method with the largest bond dimension $D=50000$, which gives reliable results with extrapolations with truncation errors in the DMRG calculation.

Figure 2: (a) Ground state energy per site $E_{0}$ for the $J_{1}$-$J_{2}$ model with $J_{2}=0.57$ and system size $14\times 14$ under OBC, as a function of the truncation errors $\varepsilon$. The data is calculated using both FAMPS and pure DMRG methods. The dashed line indicates a linear fitting, which is extrapolated to $\varepsilon=0$. Though the extrapolated results for DMRG and FAMPS are consistent, the error bar for FAMPS result is quite smaller. (b) $E_{0}$ as a function of the bond dimension $D$. For each $D$, the energy computed using FAMPS is consistently lower than that obtained with pure DMRG method. Analysis indicates that the performance of FAMPS at a given bond dimension $D$ is approximately equivalent to pure DMRG calculation with bond dimension $3D$.

To show the improvement of FAMPS over pure DMRG method, we calculate the ground state energy per site of the $J_{1}$-$J_{2}$ square lattice Heisenberg model at $J_{2}=0.57$ for system size $L=14\times 14$ under open boundary conditions (OBC). The results are presented in Fig. 2. By applying a linear extrapolation to obtain the ground state energy $E_{0}$ at the truncation errors $\varepsilon=0$, we observe that both FAMPS and DMRG methods yield consistent results. However, the FAMPS method provides a significantly smaller error bar, indicating a more reliable and accurate estimation of $E_{0}$. The FAMPS method achieves convergence more rapidly than the pure DMRG approach. Notably, even the data points corresponding to the largest truncation errors in FAMPS calculations align with the linear fitting line. In contrast, the data points with large truncation errors from the pure DMRG calculations deviate from the fitting line, indicating less robustness in approaching the ground state.

The convergence of the energy per site with bond dimension $D$ is depicted in Fig. 2 (b). The results highlight that for a given bond dimension, the energy calculated using the FAMPS method consistently falls below that calculated by the pure DMRG method. Notably, the simulation accuracy of FAMPS at a specific bond dimension $D$ is approximately equivalent to that of DMRG with bond dimension $3D$. For instance, the energy value obtained with FAMPS at $D=15000$ is comparable to those calculated with DMRG at $D=45000$.

III Results

Figure 3: (a) PVBS state with $D_{x}$ equals to $D_{y}$. (b) CVBS state with $D_{y}$ is non-zero and $D_{x}$ is vanishing.

The VBS dimer order parameter, denoted as $D_{\alpha}$ ($\alpha=x/y$), serves as a critical observable for distinguishing between the competing VBS states in the $J_{1}$-$J_{2}$ square lattice Heisenberg model Liu et al. (2018); Zhao et al. (2020). This parameter effectively captures the possible rotational symmetry-breaking, thereby providing a quantitative measure to differentiate between the CVBS and PVBS states. The definition of $D_{\alpha}$ ($\alpha=x/y$) is given by

$$D_{\alpha}=\frac{1}{N_{b}}\sum_{r}e^{iq_{\alpha}\cdot r}\hat{S}_{r}\cdot\hat{S}_{r+\alpha}$$ (5)

where $N_{b}$ is the number of bonds in the calculation, $q_{\alpha}$ is the wave vector associated with the VBS order, taking values $(\pi,0)$ or $(0,\pi)$ depending on the direction $\alpha$ of the dimer order. In the thermodynamic limit, for the CVBS state, either $D_{x}$ or $D_{y}$ will be non-zero while the other approaches zero. Conversely, for the PVBS state, both $D_{x}$ and $D_{y}$ should exhibit non-zero values, and the ratio $|\langle D_{x}\rangle|/|\langle D_{y}\rangle|$ should be unity, indicating an isotropic dimer distribution across the lattice. The illustration of the PVBS and CVBS dimer is shown in Fig. 3.

Figure 4: The bond strength $|\langle\hat{S}_{r}\cdot\hat{S}_{r+\alpha}\rangle|$ of the $12\times 12$ system at $\Delta_{y}=0.03$. The blue line shows the strong coupling bonds and the red line shows the weak coupling bonds.

Figure 5: The VBS dimer order parameter $D_{\alpha}$ versus $\Delta_{y}$ with $D_{y}$ denoted as star markers and $D_{x}$ denoted as circle markers. All data have been extrapolated to truncation errors $\varepsilon=0$. The red markers represent values extrapolated to the thermodynamic limit ($L\rightarrow\infty$). Quadratic fitting is applied to these extrapolated points versus $\Delta_{y}$. The convergence of fitting lines at $\Delta_{y}=0$ strongly suggests that there is no spontaneous rotational symmetry breaking in the VBS phase and the VBS phase is a PVBS type.

In our study, by adding an anisotropy $\Delta_{y}$ in the nearest neighboring coupling in the $y$-direction to the $J_{1}$-$J_{2}$ square lattice Heisenberg model, we can detect the possible spontaneous rotational symmetry breaking of the VBS state. As shown in Fig. 4, though the strong coupling bonds (with larger $|\langle\hat{S}_{r}\cdot\hat{S}_{r+\alpha}\rangle|$) form a plaquette structure, the rotational symmetry is broken due to the inclusion of the $\Delta_{y}$ term in Eq. (2). This rotational symmetry breaking can be detected by the VBS dimer order parameter.

In Fig. 5, we present the main results about VBS dimer order parameter $D_{\alpha}$ as a function of the anisotropy $\Delta_{y}$ for $J_{2}=0.57$ across various system sizes $L\times L$. To reduce the finite-size effect, we only use results in the central half of the studied systems. We can find that as $\Delta_{y}$ increases, $D_{x}$ ($D_{y}$) decreases (increases). For $\Delta_{y}\gtrsim 0.2$, $D_{x}$ almost approaches zero while $D_{y}$ has a finite values. This behavior is consistent with theoretical expectations of Eq. (2), where rotational symmetry is explicitly broken by including the $\Delta_{y}$ term. We also notice that as expected, when $\Delta_{y}=0$, $D_{x}$ = $D_{y}$ for given systems, because spontaneous symmetry breaking can’t occur for finite systems. To detect the possible spontaneous rotational symmetry breaking, we need to focus on how $D_{\alpha}$ behaves as $\Delta_{y}$ approaches zero.

To accurately determine the true VBS phase type at $\Delta_{y}=0$, we have to extrapolate $D_{\alpha}$ to the thermodynamic limit $L\rightarrow\infty$. In this process, the order in which the limit is taken is critical. We first approach the thermodynamic limit and then take the limit of vanishing anisotropy $\Delta_{y}$. This approach can be formally expressed as

$$D_{\alpha}=\lim_{\Delta_{y}\rightarrow 0}\lim_{L\rightarrow\infty}D_{\alpha}$$ (6)

In this study, we first perform a linear extrapolation to $L\rightarrow\infty$ for fixed $\Delta_{y}$. An example for $\Delta_{y}=0.03$ is shown in Fig. 6, and results for other $\Delta_{y}$ values are shown in Appendix. A. We subsequently apply a quadratic fitting of extrapolated $D_{\alpha}$ with $\Delta_{y}$ to obtain the spontaneous order parameters. The results are shown in Fig. 5.

Figure 6: The VBS dimer order parameter $D_{\alpha}$ versus $1/L$ for $\Delta_{y}=0.03$. The dashed line represents the linear fitting to the data points. For $D_{x}$, we use data points from $6\times 6$ to $14\times 14$, while for $D_{y}$, we exclude the data point from the $14\times 14$ system because it is not converged in the FAMPS calculation (details can be found in Appendix. B).

If the VBS phase was a CVBS type, $D_{x}$ would vanish and $D_{y}$ remain non-zero, due to the applied anisotropy $\Delta_{y}$ in the $y$-direction of the nearest neighboring coupling. But we find that the fitting curves for $D_{x}$ and $D_{y}$ exhibit remarkable convergence at $\Delta_{y}=0$, as depicted in Fig. 5. This convergence strongly suggests that the VBS phase in the $J_{1}$-$J_{2}$ square lattice Heisenberg model is a PVBS type.

One may notice that in the extrapolation process to $L\rightarrow\infty$, we utilize data points from the $14\times 14$ system size for $D_{x}$ but exclude them for $D_{y}$ (see Fig. 6). The reason for this procedure is as follows. The $D_{y}$ values obtained from the $14\times 14$ system appear unreasonably high. Moreover, the linear and quadratic extrapolations of $D_{y}$ with truncation errors in FAMPS give inconsistent values, suggesting a significant finite-bond effect for $D_{y}$ in the FAMPS calculation. However, the behavior of $D_{x}$ with truncation errors in FAMPS is quite flat and the extrapolation with truncation errors is robust. More details of the finite-bond effect are discussed in Appendix. B.

In summary, our results indicate that the VBS phase in the $J_{1}$-$J_{2}$ square lattice Heisenberg model is a PVBS type.

IV Conclusion

In this work, we employ FAMPS and DMRG methods to the $J_{1}$-$J_{2}$ square lattice Heisenberg model to investigate the nature of the VBS phase in its phase diagram. By introducing an anisotropy parameter $\Delta_{y}$ in the nearest neighboring coupling in the $y$-direction into the system, we are able to distinguish between the CVBS and PVBS types. Through careful extrapolation of truncation errors and employing reliable finite-size scaling, followed by finite $\Delta_{y}$ scaling analysis of the VBS dimer order parameter, we establish that the VBS phase in the $J_{1}$-$J_{2}$ square lattice Heisenberg model is a PVBS type, meaning the rotational symmetry is not spontaneously broken in the VBS phase.

Our results not only resolve the long-standing issue on the characterization of the VBS phase in the phase diagram of the $J_{1}$-$J_{2}$ square lattice Heisenberg model, but also demonstrate the potential of FAMPS to be applied to other complex two-dimensional systems, where large system sizes and bond dimensions are necessary to capture the intricate details of the ground state.

Appendix A Finite-size scaling

In Fig. 7, we show the extrapolation of the VBS dimer order parameter with respect to the system sizes at $\Delta_{y}=0.01,0.02,0.05$, and $0.07$.

Figure 7: VBS dimer order parameter $D_{\alpha}$ plotted as a function of system size $1/L$. The dashed line represents a linear fitting to the data points. (a) results for $\Delta_{y}=0.01$, (b) results for $\Delta_{y}=0.02$, (c) results for $\Delta_{y}=0.05$, and (d) results for $\Delta_{y}=0.07$.

Appendix B Finite-bond effect of FAMPS results for the $14\times 14$ system

To elucidate the finite-bond effect on the VBS dimer order parameter, particularly the unreliability of the results for $D_{y}$ at system size $14\times 14$, we present the extrapolation of the VBS dimer order parameter with respect to the truncation errors in FAMPS calculations for the system with size $14\times 14$ in Fig. 8. For both $D_{x}$ and $D_{y}$, linear and quadratic fittings are employed to extrapolate the truncation errors to $\varepsilon=0$. Notably, both fittings give consistent results for $D_{x}$. However, for $D_{y}$, linear and quadratic fittings give different results.

Due to the computational resource limitation, achieving the required bond dimension for the convergence of $D_{y}$ in the $14\times 14$ system size is challenging. So we exclude the data for size $14\times 14$ when performing the finite-size scalings for $D_{y}$.

Figure 8: VBS dimer order parameter $D_{\alpha}$ plotted as a function of truncation errors $\varepsilon$ for the $14\times 14$ system. The solid line indicates a linear fitting to the data points, while the dashed line represents a quadratic fitting. The blue dot-dash line represents the extrapolated value of $D_{y}$ for the $12\times 12$ system as a reference. (a) $\Delta_{y}=0.01$, (b) $\Delta_{y}=0.02$, (c)$\Delta_{y}=0.03$ and (d) $\Delta_{y}=0.05$. We can find that the extrapolation for $D_{x}$ is robust, whereas it is not reliable for $D_{y}$.