Determination of ground states of one-dimensional quantum systems using the cluster iTEBD method

TL;DR

This paper introduces a cluster iTEBD method that improves the accuracy of ground state calculations in one-dimensional quantum systems by incorporating multiple degrees of freedom, enhancing entanglement representation without increasing computational complexity.

Contribution

The authors develop a cluster version of the iTEBD algorithm that better captures quantum correlations in 1D systems, addressing accuracy issues in strongly correlated physics.

Findings

Enhanced accuracy in ground-state energy and entanglement entropy.

Successful application to spin chains with complex phase transitions.

Maintains computational efficiency comparable to standard methods.

Abstract

Within the framework of imaginary-time evolution for matrix product states, we introduce a cluster version of the infinite time-evolving block decimation algorithm for simulating quantum many-body systems, addressing the computational accuracy challenges in strongly correlated physics. By redefining the wave-function ansatz to incorporate multiple physical degrees of freedom, we enhance the representation of entanglement, thereby improving the accuracy of the ground states. Utilizing the Trotter-Suzuki decomposition and optimized truncation schemes, our method maintains roughly the same computational complexity while capturing more quantum correlations. We apply this approach to three nontrivial cases: the gapless spin-1/2 Heisenberg chain, the spin-1 anisotropic XXZD chain with a higher-order Gaussian-type phase transition, and a spin-1/2 twisted triangular prism hosting a magnetic plateau phase. Improved accuracy in physical quantities, such as magnetization, ground-state energy, and entanglement entropy, has been demonstrated. This method provides a scalable framework for studying complex quantum systems with high precision, making it suitable for situations where a pure increase in bond dimension alone cannot guarantee satisfactory results.

Equations

$$e^{-\tau H}\approx \prod _{i\in even}{e}^{-\tau h_{i,i+1}}\prod _{i\in odd}{e}^{-\tau h_{i,i+1}}.$$

$$|{\psi }_0⟩=\text{Tr}\prod _{i\in odd}{\lambda }_a^{(i)}{A}_i[a_i]{\lambda }_b^{(i+1)}{B}_{i+1}[a_{i+1}]|...a_i...⟩,$$

$$H=\sum _{\alpha }\sum _{i\in \alpha ,j\in \alpha }{h}_{i,j}^{(\alpha )}+\sum _{⟨\alpha ,\beta ⟩}\sum _{i\in \alpha ,j\in \beta }{h}_{i,j}^{(\alpha ,\beta )},$$

$$e^{-\tau H}\approx \prod _{\alpha }{e}^{-\tau h^{(\alpha )}}\prod _{\alpha \in odd}{e}^{-\tau h^{(\alpha ,\alpha +1)}}\prod _{\alpha \in even}{e}^{-\tau h^{(\alpha ,\alpha +1)}},$$

$$|{\psi }_0⟩=\text{Tr}\prod _{\alpha \in odd}{\lambda }_a^{(\alpha )}{A}_{\alpha }[m_{\alpha }]{\lambda }_b^{(\alpha +1)}{B}_{\alpha +1}[m_{\alpha +1}]|...m_{\alpha }...⟩,$$

$$A_{i,j,m_1}^{\prime }{\lambda }_{a,i}$$

$$B_{j,n,m_2}^{\prime }{\lambda }_{a,n}$$

$$\sum _{k,a_4,a_5}{R}_{i,a_{4}k}{\lambda }_{b,k}{L}_{j,a_{5}k}{G}_{a_4{a}_5,a_4^{\prime }{a}_5^{\prime }}\approx \sum _{x\le \chi }{U}_{i{a}_4^{\prime },x}{\lambda }_{b,x}^{\prime }{V}_{j{a}_5^{\prime },x},$$

$$H=$$

$$H=J_1\sum _{i;j\ne k}{\mathbf{S}}_{i,j}\cdot {\mathbf{S}}_{i,k}+J_2\sum _{i;j\ne k}{\mathbf{S}}_{i,j}\cdot {\mathbf{S}}_{i+1,k}-h\sum _{i,j}{S}_{i,j}^z,$$

$$A_{i,j,m_1}^{\prime }{\lambda }_{a,i}$$

$$B_{j,n,m_2}^{\prime }{\lambda }_{a,n}$$

Full text

I Introduction

The simulation of quantum many-body systems in strongly correlated physics poses significant computational challenges due to the exponential growth of the Hilbert space with system size. In the past decades, efficient algorithms with only algebraic complexity have been developed, such as density matrix renormalization group (DMRG) [1, 2] and time-evolving block decimation (TEBD) [3, 4]. Though designed mainly for one-dimensional quantum systems, these methods have been successfully extended to more complex geometries, such as cylinders, torus shapes, and even two-dimensional systems [5, 6]. Particularly, extending the underlying concept of matrix product state (MPS) to higher dimension [7, 8] and employing the translation symmetry, TEBD algorithms have shown great potential in the framework of imaginary-time evolution in dealing with infinite two-dimensional quantum systems, such as topological order [9, 10, 11], quantum spin liquid [12, 13], superconductivity [14, 15], and so on.

In the infinite TEBD (iTEBD) algorithm [4], the bond dimension ($\chi$) of the MPS representation of the ground state is an important hyperparameter that controls the number of variables, the upper-bound of the captured entanglement entropy, and thus the accuracy of the calculation. Therefore, in addition to designing more accurate update methods, such as canonical-form-based truncation and energy minimization variation, maintaining a larger bond dimension appears to be the only way to improve accuracy, especially when there is no suitable symmetry in the targeted systems. Therefore, it is valuable to develop a more efficient iTEBD algorithm whose performance can be improved from multiple aspects.

In this work, we address this issue by introducing another hyperparameter, the cluster size ($n$) of the MPS ansatz, to the original iTEBD algorithm. In this modified version, the ground state is represented as a special MPS form, where each local tensor supports a local physical space corresponding to $n$ particles in an ordinary MPS representation. Meanwhile, the imaginary-time evolution operator $e^{-\tau H}$ is decomposed accordingly into local terms categorized as intra-cluster and inter-cluster parts, respectively. Then, the MPS is evolved using these local terms based on the singular value decomposition (SVD) truncation strategy, assisted by entanglement spectra, until convergence is reached at an accepted precision. The algorithm has a computational cost scaling as $O(\chi^{2}d^{n+2})+O(\chi^{3}d^{n+1})+O(\chi^{3}d^{6})$, where $d$ is the dimension of the Hilbert space spanned by each particle. This complexity is the same $O(\chi^{3})$ in the large-$\chi$ limit as in the original iTEBD algorithm [3, 4], as long as $n$ is practically finite.

The validity of this cluster version of iTEBD, referred to as cluster iTEBD henceforth in this work, is demonstrated in the spin-1/2 Heisenberg chain, where better ground state energy and magnetization density are obtained, compared to the original iTEBD with the same $\chi$. We then applied this modified algorithm to two other interesting models. One is the spin-1 XXZD chain, and the other is the spin-1/2 twisted triangular prism. We verify a third-order phase transition between a Haldane phase and a symmetry-breaking large-$D$ phase in the XXZD chain, and confirm a $1/3$-plateau phase surrounded by gapless phases in the prism model. In all these three cases, purely increasing the bond dimension $\chi$ cannot guarantee satisfying results, and this is the exact situation where the cluster size $n$ can play an important complementary part.

The rest of the paper is organized as follows. In Sec. II, we briefly review the original iTEBD algorithm and introduce the modified algorithm. In Sec. III, the results of the three models mentioned above are presented in detail. In Sec. IV, we summarize and discuss some related topics.

II Method

In this section, we present the details of the algorithm proposed in this work. For this purpose, we first need to provide a brief review of the original iTEBD method, as proposed in Ref. [4].

II.1 Original iTEBD method

Employing the translation symmetry, the iTEBD method can be applied directly to the thermodynamic limit. Taking a one-dimensional translation-invariant system with only nearest neighbor interactions as an example, $H=\sum_{i}h_{i,i+1}$, for any tiny time step $\tau$, the imaginary-time evolution operator $e^{-\tau H}$ can be approximated by the Trotter-Suzuki decomposition [16],

[TABLE]

The advantage of Eq. (1) is that if the initial state $|\psi_{0}\rangle$ is represented as an MPS, described below, then $e^{-\tau H}|\psi_{0}\rangle$ can be evaluated efficiently with only local operations, and finally the ground state $|\psi_{g}\rangle$ can be obtained by applying $e^{-\tau H}$ repeatedly until some convergence is reached.

Specifically, to compute $e^{-\tau H}|\psi_{0}\rangle$, the original iTEBD method represents $|\psi_{0}\rangle$ as an MPS [17, 18],

[TABLE]

where we have assumed a two-sublattice structure, and $A_{i}$ and $B_{i}$ are local tensors supporting the physical space of the $i$-th particle. $\lambda^{(i)}_{a/b}$ is a positive vector defined on the $i$-th (odd or even) bond link, and its square mimics the entanglement spectra for the corresponding bipartition. Due to translation symmetry, the iTEBD method adopts an update strategy based on $\lambda$-assisted SVD truncation. The process is illustrated in Fig. 1. Firstly, construct a local system by grouping a local operator $e^{-\tau h_{i,i+1}}$ for the even bond (e.g., $A$-$B$ bond), two local tensors $A$ and $B$ sharing the even bond, and the corresponding odd bond vector $\lambda_{a}$ and even bond vector $\lambda_{b}$, as shown in Fig. 1(a). Then, bipartition the local system through SVD, from which the even bond vector $\lambda_{b}$ residing on the internal link of the cluster is updated as the obtained singular values. The local tensors are updated by the singular vectors appropriately, as shown in Fig. 1(b). It is worth mentioning that in the above process, the odd bond vector $\lambda_{a}$ at the edge of the cluster remains unchanged, and it will be updated similarly when evolution for the odd bond (e.g., $B$-$A$ bond) is performed, as illustrated in Fig. 1(c-d). After each SVD, to keep the bond dimension $\chi$ finite, the updated local tensors and entanglement spectra are truncated according to the singular values.

The above procedure completes an iteration of the iTEBD algorithm. This iteration is performed repeatedly until some convergence of the entanglement spectra is reached, and the converged MPS can be regarded as an estimation of the actual ground state $|\psi_{g}\rangle$. This algorithm is very efficient, and from the perspective of DMRG, the bond vector $\lambda$ is usually regarded as effective environment of the target system [18] and plays a similar role to the effective mean field in the belief propagation algorithm [19]. This idea has been extended to the two-dimensional lattice under the name simple update [6] and the mean-field version of the second renormalization group [20].

II.2 Cluster iTEBD method

In this work, we introduce another usable hyperparameter, namely cluster size $n$, to the iTEBD method, making it more flexible. Firstly, we divide the system into local clusters composed of $n$ particles, and a general local Hamiltonian can be recast as

[TABLE]

where $\langle\alpha,\beta\rangle$ means nearest neighbor cluster. The first part encompasses all interactions between particles residing within the same cluster, while the second part comprises interactions between neighboring clusters. Accordingly, the time evolution operator for a tiny Trotter step is approximated as a product of three parts,

[TABLE]

where $h^{(\alpha)}$ and $h^{(\alpha,\alpha+1)}$ contain all the intra-cluster interactions and the inter-cluster interactions, respectively.

In order to evaluate $e^{-\tau H}|\psi_{0}\rangle$, we similarly represent $|\psi_{0}\rangle$ as an MPS but with clustered physical indices, i.e.,

[TABLE]

which is the same as Eq. (2) but replacing degrees of freedom $a_{i}$ for each particle by $m_{\alpha}$ for each cluster. Suppose $d$ is the dimension of the local space for each particle, and there are $n$ particles in each cluster, then the dimension supported by $m_{\alpha}$ is $d^{n}$. In this sense, the original MPS representation expressed in Eq. (2) can be regarded as a special case, i.e., $n=1$.

As long as Eqs. (4) and (5) are established, the iTEBD method can proceed with little extra effort compared to the original version. As illustrated in Fig. 2(a), the evolution with the first part in Eq. (4) can be performed conveniently by absorbing $e^{-\tau h^{(\alpha)}}$ in the corresponding cluster directly. However, the evolution in the second part involves interactions between two neighboring clusters connected by an odd link, which requires some extra effort. Suppose $n=4$, and we need to evaluate $e^{-\tau h_{4,5}}$, where $a_{4}$ and $a_{5}$ are grouped into cluster indices $m_{1}$ and $m_{2}$, respectively. In order to save cost, as shown in Fig. 2(b), we construct a subsystem by absorbing $\lambda_{a}$ into local tensor $A^{\prime}(B^{\prime})$, and perform QR decomposition in the following

[TABLE]

where $P$ and $Q$ are column-orthogonal matrices, and $L$ and $R$ are upper-triangular matrices. Then group $L$ and $R$ as a local system and perform SVD after the time evolution, i.e.,

[TABLE]

where $e^{-\tau h_{4,5}}$ is written explicitly as its matrix representation $G$, and the index $x$ is truncated according to the generated singular value $\lambda^{\prime}_{b}$. As shown in Fig. 2(c-d), $\lambda^{\prime}_{b}$ is regarded as the updated bond vector, and local tensors are updated by multiplying the inverse bond vector $1/{\lambda_{a}}$, the column-orthogonal matrices $Q$ ($P$), and $U$ ($V$) together. One can perform the evolution with the third part involving interactions between two neighboring clusters connected by an even link, similarly.

The above procedure completes an iteration of the cluster iTEBD algorithm. Besides the bond dimension $\chi$, it provides another hyperparameter, namely the cluster size $n$, which can be used to improve the accuracy. In some sense, the MPS ansatz used in Eq. (5) can be regarded as a one-dimensional version of the projected entangled simplex states [21] proposed for higher-dimensional quantum systems, which are proved to be able to grasp more many-body entanglement than the conventional tensor network state in the context of imaginary-time evolution [22]. In the following, we will see that for a given $\chi$, within the framework of imaginary-time evolution, the MPS produced by the cluster iTEBD algorithm can capture more entanglement entropy than that obtained by the original iTEBD, and this could be the underlying reason why it can demonstrate better performance in our calculations.

III Results

In this section, the cluster iTEBD method is tested in three nontrivial models. We first validate the method through the spin-1/2 Heisenberg chain. Compared to the original iTEBD method with the same bond dimension, we obtain a smaller magnetization density and a lower ground state energy. Then, we apply it to two interesting models, where we successfully verify a third-order topological quantum phase transition and confirm a $1/3$-plateau phase, respectively.

III.1 Spin-1/2 Heisenberg Chain

The spin-1/2 antiferromagnetic Heisenberg chain has been extensively studied long before. It is well known that the exact ground state energy is $E_{\textrm{exact}}=1/4-\ln{2}=-0.443147\cdots$ from Bethe ansatz [23], and that due to strong quantum fluctuations, the ground state is magnetically disordered from the Mermin-Wagner theorem [24]. However, since the system is gapless [25], it is difficult to have an exact MPS representation of the ground state with finite bond dimension $\chi$. If SU(2) symmetry is not utilized explicitly, the MPS representation obtained from the iTEBD algorithm with a finite $\chi$ always possesses a weak but nonzero staggered magnetization. Similar phenomena have also been reported in two-dimensional systems [21, 26]. Therefore, this model provides an ideal platform to test whether the hyperparameter $n$ introduced in the cluster iTEBD method can improve the performance.

The results are displayed in Fig. 3. The smallest $\tau$ we used is less than $10^{-3}$ in our calculations. It shows that for any given cluster size $n$, the magnetization $m$ is indeed nonzero but decreases as the bond dimension $\chi$ becomes larger. Things are similar for the ground state energy $E_{g}$, i.e., for any given $n$, the relative error $\delta E_{g}\equiv\frac{E_{g}-E(\chi,n)}{E_{g}}$ decreases as $\chi$ becomes larger. More importantly, for a given $\chi$, the accuracy of the magnetization and the energy is systematically improved when $n$ is increased. This means that the cluster size $n$ indeed provides another usable control parameter to improve the performance of the iTEBD algorithm. This should be particularly useful when $\chi$ cannot be easily extended, e.g., in gapless systems where desirable $\chi$ might be quite large due to the area law of entanglement entropy [27], or two-dimensional quantum systems where large $\chi$ is difficult to handle [28].

Furthermore, in Fig. 4(a-b), we extrapolate the data of $m$ and $\delta E_{g}$ to the large-$\chi$ limit for any fixed $n$. It shows that, up to $\chi=100$ and $n=8$, the cluster iTEBD algorithm can obtain a magnetization smaller than $0.02$ and energy accuracy about $5\times 10^{-7}$, and the extrapolation for $n=8$ can give a much more accurate estimation of these quantities than that for $n=1$, i.e., the original iTEBD case. See the dashed lines in Fig. 4(a). This can be understood within the framework of imaginary-time evolution, since the cluster algorithm with the same bond dimension $\chi$ can produce an MPS that captures more quantum entanglement than the original one. The comparison of the captured entropy as a function of $n$ is demonstrated in Fig. 4(c). In practical calculations, if necessary, one can further extrapolate the obtained data $E(\chi=\infty,n)$ to the large-$n$ limit to obtain more accurate results.

III.2 Spin-1 Anisotropic XXZD Chain

We now consider the spin-1 anisotropic XXZD chain. The model has XXZ-type interactions between spin-1 operators $\boldsymbol{S}_{i}$ defined on site $i$. Moreover, the uniaxial and rhombic single-ion anisotropies are controlled by the strengths $D$ and $E$, respectively. The Hamiltonian is given by,

[TABLE]

When $E=0$, from theoretical analysis and numerical evidence [29, 30, 31, 32], it is known that the topological phase transition from the Haldane phase to a so-called large-$D$ symmetry-breaking phase can be higher-order Gaussian type, e.g., third-order for $J_{z}=1$ at $D_{c}\approx 0.9687$ and fifth-order for $J_{z}=0.5$ at $D_{c}\approx 0.6197$ [32]. Meanwhile, if we fix $D=D_{c}$ and vary $E$, there is an ordinary second-order phase transition between the two Néel phases breaking $Z_{2}$ symmetry, e.g., at $E=0$ when $J_{z}=1$. A sketch of the phase diagram for $J_{z}=1$ is illustrated in Fig. 5. This model has drawn attention recently, partly due to its relevance to a recent ultra-cold atom realization [33] and inelastic neutron scattering experiments on real materials [34].

In this work, we set $J_{z}=1$ for simplicity, and focus on the two phase transitions, namely the third-order topological phase transition and the second-order symmetry-breaking one, as indicated by blue and orange dashed lines in Fig. 5, respectively. The third-order transition is not easily accessible. The previous calculations utilize complex quantities, such as fidelity [29], level crossing point [31], or cross derivatives [32], to explore the signal, as the conventional second-order derivatives of the energy cannot yield a clear peak [30]. In this work, we expect a better identification of peaks in the usual second-order derivatives by increasing the hyperparameter $n$.

The results are summarized in Fig. 6 and Fig. 7. In Fig. 6, setting $E=0$, we plot the ground state energy $E_{g}$ as a function of $D$ as well as its derivative $\frac{\partial E_{g}}{\partial D}$ and second derivative $\frac{\partial^{2}E_{g}}{\partial D^{2}}$. It shows that, indeed, there are no obvious singularities in the curves; instead, the $\frac{\partial^{2}E_{g}}{\partial D^{2}}$ develops a broad peak that becomes more and more significant as either $n$ or $\chi$ increases. The location of the peak for $(n=6,\chi=100)$ is about $D\approx 0.964$, providing a consistent estimation of the third-order phase transition with previous studies [32]. On the contrary, as shown in Fig. 7, by fixing $D=D_{c}$ obtained for each ($n,\chi$) setup and plotting the derivatives with respect to the rhombic anisotropy $E$, we can see clear singularities of $\frac{\partial^{2}E_{g}}{\partial E^{2}}$ at $E=0$ for a series of $(n,\chi)$ combinations. This is consistent with the orange dashed line in Fig. 5 indicating the second-order phase transition.

We can see from this calculation that for the ordinary continuous phase transition with symmetry breaking which can be determined by the original iTEBD algorithm, the cluster iTEBD can provide more accurate critical information, such as the critical value $D_{c}$ as shown in Fig. 6(c) and indicated in the legends of Fig. 7. While for the Gaussian-type transition, where there is no singularity in the second-order derivatives, the iTEBD algorithm with $n=1$ cannot determine the phase transition easily, and the cluster iTEBD with larger $n$ can display peaks more clearly.

III.3 Spin-1/2 Twisted Triangular Prism

The third model we studied is a spin-1/2 twisted triangular prism, which can be effectively depicted in Fig. 8. The triangular prism is composed of three spin chains (horizontal lines), and the three spins on the same cross-section construct a triangle. The spins on the same triangle are assigned different colors, and by twist, we mean that each spin interacts with only those of different colors in the same or nearest neighbor triangles. The Hamiltonian can be effectively expressed as

[TABLE]

where $\boldsymbol{S}_{i,j}$ is the spin-1/2 operator defined on the $i$-th triangular, and $j=1,2,3$ denotes its color index. $J_{1}$ means interaction strength in the same triangle (denoted by black in Fig. 8), $J_{2}$ denotes interaction strength between nearest neighbor triangles (denoted by blue in Fig. 8), and $h$ represents a uniform external magnetic field in the $z$-direction. Obviously, if the system is considered a one-chain system, there must exist interactions between at least third-nearest neighbors, which is not easy for the conventional iTEBD to handle.

In this work, we set the energy scale as $J_{1}=1$ for simplicity. Due to the strong frustration, the model hosts a rich phase diagram in the range $0<J_{2}\lesssim 1.22$ [35, 36]. Specifically, when turning on $h$ within this parameter range, as determined by DMRG and exact diagonalization on finite systems, the model was found to host a $1/3$-plateau state surrounded by gapless magnetically ordered phases [35]. The model is believed to be relevant to some real compounds, such as $[(\textrm{CuCl}_{2}\textrm{tachH})_{3}\textrm{Cl}]\textrm{Cl}_{2}$ [37], and it is argued to apply to other frustrated spin tubes as well.

Using the cluster iTEBD algorithm, the MPS ansatz in Eq. (5) suitable for this system can be conveniently represented by grouping the spin degrees of freedom in each triangle as a single composite particle, which spans a local Hilbert space with dimension $2^{3}=8$. Then, the intra-triangle and inter-triangle interactions are reduced to just on-site and nearest-neighbor ones, respectively. In this work, we use the minimal $n=3$ for this model in all the calculations.

We focus on two specific parameters, namely $J_{2}=0.1$ and $J_{2}=1.2$, which are close to the two edges of the interesting range mentioned above. The results are illustrated in Fig. 9 and Fig. 10. For $J_{2}=0.1$, we calculate the ground state energy $E_{g}$, its derivative $\frac{\partial E_{g}}{\partial h}$, as well as magnetization $m_{z}$, as functions of magnetic field $h$. It shows that the derivative is a constant over a wide range about $0.195\lesssim h\lesssim 1.51$, and this is consistent with the fact that the magnetization displays a $\frac{m}{1/2}\approx\frac{1}{3}$ plateau at the same range. The plateau in $\frac{\partial E_{g}}{\partial h}$ at $h\gtrsim 1.8$ corresponds to the fully polarized state residing at the magnetic plateau with $m_{z}=0.5$. For $J_{2}=1.2$, the derivative $\frac{\partial E_{g}}{\partial h}$ and $m_{z}$ also show plateau behavior but at the range of $1.55\lesssim h\lesssim 2.55$. Furthermore, to reveal the nature of the $1/3$-plateau state, in Fig. 10(c), we plot the magnetization density for each spin separately in the same triangle. We find that the spins denoted by black and red ($\boldsymbol{S}_{i,2}$ and $\boldsymbol{S}_{i,3}$) tend to be parallel to each other, while spins colored green ($\boldsymbol{S}_{i,1}$) are always anti-parallel to them. This should correspond to the so-called UUD state found in Ref. [36], and we do not find a plateau state with uniform magnetization claimed there.

In fact, in Fig. 9, besides the $1/3$-plateau and saturated plateau, there also seems to exist a $1/9$-plateau phase at the range of $0.13\lesssim h\lesssim 0.14$. Similar things happen for $J_{2}=1.2$ at about $0.8\lesssim h\lesssim 0.9$. To verify this possibility, we extend our calculation to $\chi=150$, and the result is summarized in Fig. 11(a) for $J_{2}=1.2$. It shows clearly that as $\chi$ becomes larger, the width $w$ of the plateau becomes narrower and narrower. This is confirmed further by the entanglement entropy extracted from the MPS representation of the ground state, as shown in Fig. 11(b): the entanglement entropy is strongly suppressed in the plateau phase, and the range of the suppressed region becomes smaller as $\chi$ becomes larger. What’s more, in the inset of Fig. 11(a), we find that the data of $w$ can be well fitted by linear functions of $1/\chi$, and the extrapolation to the large-$\chi$ limit gives $w\simeq 0$. Therefore, the obtained $1/9$-plateau is probably a metastable state arising from the finite-$\chi$ effect and should be absent in the true ground state.

IV Summary and Discussion

In this study, we have proposed a cluster version of the iTEBD algorithm to accurately determine the ground state of a one-dimensional quantum lattice system. Similar to the original one that has been widely used, the cluster algorithm involves designing a clustered MPS ansatz, decomposing the imaginary-time evolution operator as intra- and inter-cluster parts accordingly, and updating the MPS by entanglement-spectra-assisted SVD truncation. The cluster size $n$, as well as the bond dimension $\chi$, provides useful control parameters for improving the performance. Its validity has been demonstrated in the spin-1/2 antiferromagnetic Heisenberg chain, where the cluster iTEBD can provide a lower ground state energy and smaller magnetization compared with the original iTEBD with the same bond dimension. After that, we first apply it to the spin-1 anisotropic XXZD chain. Through the second derivatives $\frac{\partial^{2}E_{g}}{\partial D^{2}}$, we verify a third-order Gaussian-type phase transition between the Haldane phase and the large-$D$ phase, as well as a second-order phase transition between two Néel phases. Then, we apply it to the spin-1/2 twisted triangular prism, and identify a $1/3$-plateau phase and metastable $1/9$-plateau phase through magnetization and entanglement entropy calculations.

Given cluster size $n$, bond dimension $\chi$, and degrees of freedom $d$ of a single particle, the cluster iTEBD algorithm has computational cost scaling as $O\left(\chi^{2}d^{2n}\right)+O\left(\chi^{3}d^{n+1}\right)+O\left(\chi^{3}d^{6}\right)$, if neither symmetry nor partial SVD has been considered. The first part comes from intra-cluster evolution and can be reduced to $O\left(\chi^{2}d^{n+2}\right)$ by further decomposing the $e^{-\tau h^{(\alpha)}}$ bond by bond separately, instead of regarding it as a whole. The remaining two parts come from inter-cluster evolution involving Eq. (6) and Eq. (7), respectively. When $n$ is finite, in the large-$\chi$ limit, the cluster algorithm has the same order of $\chi$ as the original iTEBD, namely $\chi^{3}$.

Capturing the entanglement entropy more easily and being able to systematically improve performance by enlarging $n$, the cluster iTEBD is advantageous in at least two situations. One is the situation where purely increasing $\chi$ is difficult to achieve satisfactory results, such as fast convergence of physical quantities in gapless systems like the spin-1/2 Heisenberg chain discussed in Sec. III.1, and pronounced peaks in derivatives like the spin-1 XXZD chain studied in Sec. III.2. The other is the system hosting interactions with longer ranges than nearest neighbors, e.g., the spin-1/2 twisted triangular prism described in Sec. III.3. When dealing with longer-range interactions, such as $h_{2,6}$ for $n=4$ ansatz depicted in Fig. 2(a), to perform the time evolution of $e^{-\tau h_{26}}$, Eq. (6) can be modified similarly by separating $a_{2}$ and $a_{6}$ from other indices, i.e.,

[TABLE]

We note that the form of Eq. (7) is unchanged, but the indices $a_{4}$ and $a_{5}$ are now replaced by $a_{2}$ and $a_{6}$, respectively. Then, the subsequent procedures remain the same as before. In fact, it is easy to show that as long as the interaction range is no larger than $n$, the cluster iTEBD algorithm with cluster size $n$ involves only the nearest neighbor interactions.

The cluster iTEBD algorithm shares some resemblance with the projected entanglement simplex state method [21] and the so-called regularized time evolution scheme [26], in the sense that in these work, multiple degrees of freedom are also grouped together in either wave function ansatz or time evolution operator. It is also worth mentioning that, in fact, it has little to do with previous studies of cluster update, which updates a two-dimensional state in the imaginary-time evolution by regarding a small cluster nearby as the effective environment of a target system [38]. In general, there are no grouped degrees of freedom in the cluster update algorithm.

In this work, the maximal kept bond dimension $\chi$ is no more than 150, because it is already sufficient for our specific purpose. However, it is worth mentioning that since the cluster iTEBD method can produce an MPS with more entanglement than the original one (though the upper bound $\log\chi$ remains unchanged), as shown in Fig. 4(c), it is necessary to use larger $\chi$ to obtain the same truncation error as the original iTEBD. In other words, the cluster iTEBD algorithm provides a better estimation of how large $\chi$ should be used for a given system, and it makes no sense to use a large $n$ but with a small $\chi$. In practical calculations, one should make a balance between $n$ (entanglement) and $\chi$ (truncation error), which is similar to the case of two-dimensional quantum systems [28]. In any case, keeping a sufficient $\chi$ is always necessary, as alerted in Fig. 11 for the $1/9$ magnetic plateau in the prism model.

Acknowledgments

We are grateful to Rui-Zhen Huang, Hongyu Chen, and Rong Yu for their helpful discussions. This work was supported by the National R&D Program of China (Grants No. 2023YFA1406500, 2024YFA1408604, 2022YFA1403601, 2024YFA1410500), the National Natural Science Foundation of China (Grants No. 12274458, No. 12322402, No. 12274206), the Innovation Program for Quantum Science and Technology (Grant No. 2021ZD0302800), the Natural Science Foundation of Jiangsu Province (No. BK20233001), and the Xiaomi foundation.