1 Introduction↩︎

SU(\(N\))-symmetric quantum many-body systems are one of the central fascinating research areas in modern condensed matter physics. They offer rich phase diagrams and exotic quantum phases that extend far beyond the conventional SU(2) spin systems. The experimental realization of high spin symmetries has been achieved on several platforms, most notably with cold gases of alkaline-earth(-like) atoms, where nuclear spin degrees of freedom can be decoupled from electronic motion, effectively realizing SU(\(N\)) symmetric interactions with \(N\) up to 10 [1], [2], as well as, more recently, with ultracold molecules reaching even higher values of \(N\) up to 36 [3]. Additionally, transition metal compounds with orbital degeneracies naturally exhibit SU(4) symmetry when spin and orbital degrees of freedom are treated on an equal footing, as observed in certain vanadium and titanium oxides [4]–[6]. From a theoretical standpoint, SU(\(N\))-symmetric models provide fundamental insights into quantum criticality, entanglement properties, and the interplay between symmetry and topology in quantum many-body systems [7], [8]. The computational study of these models, however, presents significant challenges due to the exponential growth of the Hilbert space dimension with \(N\), necessitating the development of advanced numerical techniques that can efficiently capture the entanglement structure of quantum states.

Tensor network methods [9]–[11] have shown great promise in the study of strongly correlated quantum systems, including SU(\(N\))-symmetric models. These theoretical approaches approximate the ground-state wave function or thermal density matrix as a contraction of a large number of small-rank tensors—the tensor network. This network of tensors is chosen in a form that faithfully represents the structure of correlation and entanglement of the ground state. For one-dimensional (1d) quantum lattice systems, the most widely used tensor network ansatz is the matrix product states (MPS), which is also the basis of the celebrated density matrix renormalization group (DMRG) algorithm [12], [13]. The direct generalization of MPS to two-dimensional (2d) systems in the thermodynamic limit is the infinite projected entangled pair states (iPEPS) [14]–[17] (for practical introductions, see also Refs. [18], [19]). The iPEPS approaches have been successfully employed to study highly correlated 2d problems: the Heisenberg model on the kagome lattice [20], the doped Hubbard model on the square lattice [21], as well as topological orders [22]–[24]. The 2d iPEPS algorithms were extensively applied to SU(\(N\))-symmetric Heisenberg and Hubbard models on various lattices, including square [25]–[28], honeycomb [29]–[32], triangular [25], and kagome [33], [34] geometries.

The iPEPS tensor networks have recently been further extended to three-dimensional (3d) quantum systems [35]–[38]. The main difficulty of the 3d iPEPS approaches is the calculation of observables. This cannot be done exactly and should be completed using an approximate method. Several approximate solutions were proposed in the literature: first, the Refs. [36], [39], [40] suggested using Simple Update mean-field environments to calculate observables. This method can be accurate enough only for the gapped models. Second, Ref. [35] proposed either employing cluster environments or lower dimensional tensor-network methods, such as 2d boundary iPEPS and corner transfer matrix renormalization group (CTMRG), to calculate observables. In our recent work [41], we have further generalized the boundary iPEPS + CTMRG approach of Ref. [35] by employing the “single-layer” mapping for both boundary PEPS and CTMRG [20], [42]–[46]. This allowed us to significantly reduce the computational cost of the observables calculation for 3d iPEPS and reach bond dimensions as high as \(D=6\). Other relevant 3d tensor-network studies include Refs. [47]–[52].

In this paper, we apply tensor-network methods to the SU(\(4\))-symmetric Heisenberg model on the cubic lattice. We first introduce the model and discuss the results obtained previously in its 1d and 2d limits. We analyze the highly-anisotropic limits and further elaborate on the appearance of previously reported dimerized phases [26]. Finally, we analyze the three-dimensional cases, both anisotropic and fully isotropic, using the 3d iPEPS approach developed earlier in Refs. [35], [41]. Our previous work [41] focused on a rather simple 3d lattice problem: the SU(2)-symmetric Heisenberg model on the cubic lattice with \(1\times1\times1\) unit cell. To study the more complex SU(4)-symmetric Heisenberg model, we generalize the proposed method to larger unit cells. We describe the necessary generalizations and apply these to the system under study.

2 Model↩︎

The main focus of our current study is the SU(4)-symmetric Heisenberg model on the simple cubic lattice. The internal degrees of freedom on the lattice sites (sometimes also called colors) transform in the fundamental representations of the symmetry group (the local Hilbert space on each lattice site has dimension \(4\)). The Hamiltonian of this model has the following form: \[\label{eq:model} H = \sum_{\langle ij \rangle_{\gamma}} J_{\gamma} P_{i,j},\tag{1}\] where \(\langle ij \rangle_{\gamma}\) denotes summation over the nearest-neighbor pairs of sites along the crystallographic axis \(\gamma = \{x, y, z\}\) and \(P_{ij}\) is the permutation operator, which acts on the two local basis states on sites \(i\) and \(j\) as follows: \(P_{ij} |\alpha_i \beta_j\rangle = |\beta_i \alpha_j \rangle\). Note that the coupling constants \(J_{\gamma}\) can be different; however, if not stated otherwise, we fix \(J = J_{x} = J_{y} \neq J_{z}\). The case in which all couplings are equal is referred to as the 3d isotropic point.

In the limit \(J_{z} \to 0\), the layers in \(xy\)-planes decouple from one another. It is known that in this limit, the ground state of the SU(4) Heisenberg model on the square lattice spontaneously breaks both translational and cyclic rotational symmetry, as well as continuous SU(4) spin symmetry. The many-body ground state in this limit is characterized by the \(4 \times 2\) unit cell, with different colors forming pairwise dimers [26].

In the opposite limit \(J_{z} \to \infty\), the chains of sites in the \(z\) direction also decouple from each other and are described separately by the SU(4)-symmetric Heisenberg models on the chain. It is known that the one-dimensional SU(4) Heisenberg chain (completely decoupled from the other neighbors) does not break any symmetries; the system resides in the critical phase, with its low-energy theory described by the Wess-Zumino-Witten (WZW) conformal field theory [53], [54].

Our main claim in this study is that the dimerized color-ordered phase, which is known from the isotropic 2d limit, persists throughout the entire phase diagram of the model for arbitrary finite \(J_{z}/J\) and is continuously connected to the limit of decoupled WZW models, as illustrated in Fig. 1. In particular, we argue that the coupling of SU(4) chains in the transverse directions induces dimerization, which leads to the known ground-state configuration of Néel-ordered color dimers.

3 Methods↩︎

We employ several tensor-network methods to investigate the model under study. The main reason is that the SU(4)-symmetric Heisenberg models are capable of hosting various highly entangled phases, such as valence-bond solids or spin liquids [26], [30]. To properly describe these strongly correlated regimes, it is necessary to work within a theoretical framework capable of representing highly entangled states. Among such methods, the most widely used are either tensor network or variational Monte Carlo approaches (including the novel neural quantum state algorithms). In this respect, tensor networks have a slight advantage, as they can be directly applied in the thermodynamic limit. Note that, below, we generally employ real-valued tensor networks. For verification purposes, we also developed other alternative optimization schemes with complex-valued tensor networks, which resulted in the same observations. We also do not impose any additional symmetries on our tensor networks.

Our theoretical analysis consists of three parts. First, we briefly return to the SU(4)-symmetric Heisenberg model on the square lattice, which was already studied in Ref. [26]. The difference from the previous study is that here we also consider the anisotropic cases with \(J_{x} \neq J_{y}\). We employ 2d iPEPS to obtain new results in this regime. Second, we consider another limit, where \(J_{z} \gg J_{x}=J_{y}\) and work in the framework of the decoupled MPS (decoupled chains) ansatz to study this regime with the MPS plus mean-field theory (MPS\(\,+\,\)MF) approach [55]–[57]. Third, we focus on the genuine 3d isotropic point with \(J_{x}=J_{y}=J_{z}\) and employ 3d iPEPS. These three cases are described in detail below.

3.1 2d iPEPS↩︎

We begin our theoretical analysis from the SU(4)-symmetric Heisenberg model on the square lattice with anisotropic couplings (\(J_{x} \neq J_{y}\)). For this purpose, we employ the standard iPEPS ansatz, as shown in Fig. 2 (b). We take the elementary unit cells of size \(4 \times 2\) and \(4 \times 4\). The iPEPS is initialized with a random product state and then optimized using the Simple Update method with the gradually decreased time step [18], [19], [58]. After the optimization, we compute energies for different lattice bonds using the double-layer CTMRG method [59], [60].

3.2 MPS + MF↩︎

We are also interested in the limit where the coupling along one crystallographic axis is much stronger than the couplings in the transverse plane. In this limit, we can view the system as a set of nearly decoupled, weakly interacting chains. We expect that the correlations in these chains are strong, as is natural for one-dimensional systems, while the entanglement between different chains should be weak. In the extreme case, we can consider that the interaction between different chains can be included by the mean-field treatment, while the entanglement along the chain should be fully described with a non-perturbative method. In this limit, we propose to use the ansatz consisting of separate matrix product states (MPS) along different chains, as shown in Fig. 2 (c). Such a combination of MPS along the chain and mean-field theory treatment of the interaction between the chains is known as the MPS + MF framework [55]–[57]. To optimize such an ansatz, we employ infinite time-evolving block decimation (iTEBD), where, after the application of imaginary time Trotterized evolution gates between different chains, we truncate the new bond back to the bond dimension \(D=1\). Note that, due to the fact that the chains are decoupled, we can compute all the observables exactly without CTMRG. During the optimization, we gradually increase the bond dimensions of separate MPS.

3.3 3d iPEPS↩︎

In the proximity of the isotropic point of the model 1 , we cannot expect decoupled chains or 2d iPEPS wave functions to be suitable ansatzes, as generally, we can expect rather complex 3d correlations. To include such a possibility, we employ 3d iPEPS wave functions, as shown in Fig. 2 (a). The 3d iPEPS wave function is optimized using imaginary time evolution with the Simple Update method [18], [19], [58], as described in Refs. [35], [36], [39], [40].

As an initial approximation, we compute observables with the Simple Update environments (see Ref. [36]). This calculation is approximate, but it allows one to estimate the phase of the wave function (perhaps, except for the states with topological order, which is not the case for the SU(4)-symmetric Heisenberg model). In addition, it allows one to pinpoint the necessary unit cell and iPEPS bond dimensions \(D\). After this initial step, the observables of iPEPS are computed with a more accurate CTMRG method, which is described below.

3.4 Boundary iPEPS and single-layer CTMRG↩︎

To compute operator expectation values, it is necessary to contract an infinite double-layer 3d tensor network, as shown in Figs. 3 (a) and (b).

Such tensor network contraction is performed approximately in two steps. First, the tensor network parts from the top and from the bottom, which sandwich a chosen horizontal layer, are contracted. Let us first assume that the unit cell has the lowest (i.e., unit) periodicity in the vertical direction. To contract the bottom and top parts of the tensor network, we should note, first, that a horizontal layer of the 3d tensor network can be viewed as a matrix formed by its top to bottom indices. Then, the contraction is equivalent to finding the leading left and right eigenvectors of this matrix. Due to the inherent lattice structure of the matrix, we can assume that the leading eigenvectors can be represented as 2d iPEPS (boundary iPEPS or bPEPS) with bond dimension \(\chi_{b}\), as shown in Fig. 3 (c). To find these leading eigenvectors, we employ the power method, as described in detail in our previous work [41]. We apply every double layer to the boundary iPEPS in two stages, as discussed in Ref. [41]. After each stage, we truncate the bond dimensions of the bPEPS indices back to their original values using superorthogonal canonical form [61], [62]. We can now weaken our initial assumption of unit periodicity in the \(z\)-direction and assume that there are \(k\) different layers stacked vertically and repeating periodically. In this case, the contraction of the tensor network from the top and bottom is equivalent to finding the leading eigenvectors of the matrix obtained as multiplication of the \(k\) layer transfer matrices. To obtain these eigenvectors, we again assume that the eigenvectors can be approximated as 2d iPEPS and then apply these \(k\) layer transfer matrices to these 2d iPEPS sequentially (similar to how it is done in 2d case with the boundary MPS approaches [63]). After application of every layer, we perform truncation of bPEPS bond dimensions.

After we determined the bPEPS tensor networks, we contract the residual four layer 2d tensor networks shown in Fig. 4. To this end, we first map the resulting tensor network onto the single-layer network, as shown in Fig. 4 and discussed in detail in Ref. [41]. The main difference in this mapping compared to Ref. [41] is that the original unit cell in the \(xy\) plane can be nontrivial (e.g., of size \(m\times n\) with \(m,n\in \mathbb{N}\)). Then, the mapping results in the single-layer tensor network with the effective unit cell \(4m \times 4n\). This new single-layer tensor network is then contracted with the CTMRG algorithm and auxiliary bond dimension \(\chi_{c}\) [59], [60], [64], [65]. Note that one may sometimes encounter the ill-convergence of this CTMRG by starting from random values. To mitigate such potential issues with convergence, we use two different approaches: (i) physical initialization of CTMRG tensors followed by a couple of layers of bMPS based contraction prior to full CTMRG; (ii) CTMRG with larger environments for the definition of projectors (in particular, if for computing projectors one usually contracts \(4 \times 2\) blocks formed by the tensor network and CTM tensors [29], then in our case, we use \(4 \times l\) blocks, where \(l>2\)).

In the discussion above, we generally assumed that the tensor network is contracted from the top and bottom of a chosen layer with bPEPS. This contraction allows us to compute operator expectation values for the operators inside this layer. Hence, the CTMRG contraction should be repeated for all the \(k\) layers. Still, the described construction cannot be used to compute the energy expectation values for operators along vertical \(z\)-bonds. To deal with such operators, we rotate the tensor network and repeat the calculations above with the rotated network, where the new horizontal layer corresponds, e.g., to the original \(xz\) planes.

Finally, there are additional self-consistency checks, which enable verification of our results. First of all, the local reduced density matrices are hermitian, but the contraction method described above does not guarantee this hermiticity, since the top and bottom bPEPS are generally different, and the single-layer CTMRG also does not use hermiticity. Hence, the hermiticity of the density matrix is restored only after sufficient convergence of bPEPS and CTMRG—this is one of the nontrivial checks of the validity of the algorithm. Second, as mentioned above, certain observables require additional tensor network rotations to be found, while other observables, such as the onsite reduced density matrices, can be determined by using both the original and rotated tensor network contractions. Therefore, we compute these onsite reduced density matrices by utilizing both mentioned contractions and compare the results to confirm the algorithm convergence and consistency.

4 Results↩︎

4.1 2d limit↩︎

We first turn our focus to the systems in the limit of decoupled planes, i.e., square lattices with anisotropic in-plane couplings. In 2d case, the full machinery of 2d iPEPS is available to us.

In particular, we analyze the 2d limit by taking \(J_{z} = 0\) and \(J_{x} \neq J_{y}\). Without loss of generality, below we assume that \(J_{x} \geq J_{y}\) and \(J_{x}=1\). In the isotropic case, it is known that the model is in the color-ordered dimerized phase [26], which spontaneously breaks lattice translational, lattice rotational and continuous spin-rotational SU(4) symmetries. In this subsection, we study the anisotropic case with \(J_{y} \in [0,1]\) by applying the iPEPS algorithm and show that the dimerized phase at \(J_{x} = J_{y}\) is continuously connected with the limit \(J_{y} \to 0\), which is described by the decoupled translationally invariant WZW models, with translational invariance [and SU(4) symmetry] along the chain gradually restored.

In particular, we take the 2d iPEPS wave functions with unit cells of size \(4 \times 2\) and bond dimensions \(D_x=8\) and \(D_y=6\) along the \(x\) and \(y\) directions, respectively, as the ground-state ansatzes. We employ the CTMRG approach with the auxiliary dimension \(\chi = 80\) for the computation of two-site correlation functions. In all cases, we observe the same structure of the dimerized phase as in the isotropic limit.

In Fig. 5, we show the energy expectation values along the \(x\)-direction (along the dimer orientation) for both strongly-correlated and weakly-correlated bonds. Note that here and below, the bond-specific energies \(E_{ij,\gamma}=\langle H_{ij,\gamma}\rangle\) are additionally averaged over all similar bonds in the given unit cell. We can conclude that the breaking of translational symmetry, as well as the continuous SU(4) symmetry, is a consequence of the addition of a (small) coupling between different chains. Moreover, the greater the interchain coupling \(J_y\), the more pronounced the breaking of translational symmetry along the chain becomes. It means that the coupling between different chains is the reason for the dimerized color-ordered phase, which remains stable including the 2d isotropic point.

4.2 1d limit↩︎

We proceed to study the limit \(J_{x} = J_{y} \to 0\) while keeping the coupling between chains along the \(z\)-direction fixed (\(J_z=1\)). The decoupled chains are described by the WZW models and are translationally invariant. As discussed in the previous subsection, the addition of a small coupling between the chains may naturally lead to dimerization. In particular, in the mean-field description, the interaction with the four nearest-neighbor chains is simply the same as the interaction with the two chains with the twice larger couplings. Hence, we can expect that the behavior of coupled chains can be qualitatively similar on both highly anisotropic square and cubic lattices (with the difference that the cubic lattice may potentially have more complex unit cells).

We study the behavior in this limit within the decoupled chain MPS + MF ansatz. Our results for the energy dependence on different lattice bonds are shown in Fig. 6 on the cubic lattice with unit cell \(4 \times 4 \times 4\). It is clear that, in both cases, dimerization develops as a consequence of interchain coupling. The results are obtained for the MPS bond dimensions up to \(D_z=70\).

We should also emphasize several details here. First, at small interchain coupling (\(J_{x,y}\lesssim0.1J_z\)), the energy per bond smoothly connects to the WZW model with a uniform distribution of energies along the chain with \(E_z\approx-0.825\) [54]. Second, the expectation values of the exchange operators \(P_{x,y}=E_{x,y}/J\) between different chains (which are computed in a mean-field-like fashion, since the neighboring chains are not entangled) decrease from \(0.25\) to \(0.1\). The limit with \(P=0.25\) corresponds to the probability that the sites in two different chains have the same color when the colors are distributed homogeneously along the chain. The decrease of this expectation value to \(0.1\), even at the level of mean-field theory, points to dimerization in the chains. Finally, we observe that for the sites on the same chain, the expectation value of the operator \(P\) differs between even and odd bonds. This also points toward the dimerization of the chains.

4.3 Isotropic point in 3d↩︎

From the two previous subsections, it is clear that dimerized color-ordered states appear in the cubic lattice both in the limit of large \(J_{z}\) and in the opposite limit \(J_{z} \to 0\). The natural question is whether these two regimes are adiabatically connected to each other through the isotropic point. In this subsection, we present the evidence that this is, indeed, the case. For this purpose, we study the system at the isotropic point within the 3d iPEPS approach.

We first optimize the iPEPS wave function with the Simple Update on the \(4 \times 4 \times 4\) lattice and with bond dimensions up to \(D=8\). To approximately pinpoint the phase, we calculate the observables using the Simple Update environments. Starting from \(D=4\), the wave function converges to the dimerized phase with \(2 \times 2 \times 4\) unit cell. After this preliminary estimation of the necessary unit cells and bond dimensions, we focus on the \(2\times2\times4\) unit cell with the bond dimension \(D=5\). The bond dimension \(D=5\) was chosen for two reasons: first, Simple Update calculations with \(D=4\) showed less pronounced convergence to the dimerized state (sometimes converging to product-like states), while for bond dimensions \(D \geq 6\), the necessary computational cost becomes too large. The results below are shown for the \(D=5\) wave function.

Next, we employ the boundary iPEPS + CTMRG approach to calculate observables of the optimized 3d iPEPS wave function. In particular, we calculate the local onsite reduced density matrices for all sites inside the unit cell, as well as the energies on the bonds. It is necessary to converge boundary iPEPS in \(z\) and \(y\) directions to calculate energies on all bonds. For \(z\) direction, the maximal boundary iPEPS bond dimension is fixed to \(\chi_{b} = 14\), while for \(y\) direction, the maximal boundary iPEPS bond dimension is fixed to \(\chi_{b}=17\). For these bond dimensions, the smallest boundary iPEPS bond spectra \(\lambda_{b}\) are of the order \(10^{-5}\), which is used as a criterion of bPEPS convergence. For converged bPEPS, we finally calculate the observables using the single-layer CTMRG, with the CTMRG bond dimensions up to \(\chi_{c}=160\). The corresponding bond energies converge reaching accuracy in the range \(10^{-5} - 10^{-4}\), while the local reduced density matrix elements converge with accuracy below \(10^{-5}\). Note that the local reduced density matrices also converge with respect to hermiticity up to errors of the order of \(10^{-6}\), and the difference between local reduced density matrix elements computed with the original and rotated bPEPS is of the order of \(10^{-5}\).

We observe that the bonds can be divided into three groups according to their respective energies: the first group consists of \(x\) and \(y\) bonds; the second group includes energies on \(z\) bonds starting from even \(z\)-coordinates; the third group includes the rest of \(z\) bonds. The average energy of bonds in the first group is equal to \(E_{x} \approx E_{y} = -0.1823\). The energies of bonds in the second group have the average energy \(E_{z, even} = -0.8668\), while the energies in the third group are equal to \(E_{z, odd} = - 0.1908\). Note that, as various rotational and translational symmetries are not enforced during the optimization and Simple Update iPEPS optimization being not very accurate, there were some residual discrepancies between energies on bonds of the same group, especially for the bonds in the first and third groups.

In Fig. 7, we show the onsite concentrations of different colors on different sites of the unit cell and also show the bonds with the largest energies (these we call above as the bonds of the second kind). The concentrations of colors on the site approximately follow the proportions \(0.447:0.442:0.057:0.054\). The slight difference between the first and second concentrations is likely due to the limitations of the Simple Update at the bond dimension \(D=5\). Still, we can conclude that the concentrations of the first two and the last two colors are approximately equal on all the sites. In Fig. 7, we can clearly see that the colors form either “red-green” or “yellow-blue” dimers with large negative energies.

5 Conclusions and Outlook↩︎

We studied many-body ground states of the SU(\(4\))-symmetric Heisenberg model on the simple cubic lattice using different tensor network approaches, in particular, by employing three-dimensional iPEPS variational wave functions. For this wave function ansatz, we extended our previous method of calculating observables to handle large unit cells. Our numerical results point towards the persistence of the dimerized color-ordered phase known from the 2d limit [26] in both anisotropic and completely isotropic 3d regimes. This, in turn, is a sign that the three-dimensional SU(\(N\))-symmetric lattice systems remain very promising as hosts for exotic quantum phases beyond the mean-field theory results [66].

There are several possible additional research directions. One of them is the extension of the method to different lattice geometries. The simplest lattices to extend our calculations are hyperhoneycomb and diamond lattices, as the tensor networks on these lattices can be easily mapped onto the simple cubic lattice after the combination of several tensors into one. It is also possible to develop algorithmic extensions to pyrochlore and octahedron lattice geometries (3d Lieb lattice) after switching from iPEPS to the infinite projected simplex states (iPESS) ansatz [67]. Some other cases, such as the face centered cubic lattice, may also be within reach after the inclusion of longer range interactions. In this respect, the most straightforward direction to proceed is to study the Heisenberg model on the pyrochlore lattice [68] and some additional models on the hyperhoneycomb lattice. Note that it was previously proposed that some SU(\(4\))-symmetric models may be realized on the hyperhoneycomb lattices [69].

Besides potential extensions to other lattice geometries, it is also possible to further develop the algorithm on the cubic lattice. One potential extension is the implementation of the Full Update optimization scheme [70], [71] for the 3d iPEPS wave functions. The cost of such optimization is higher but technically manageable. In conditions where we are restricted to small bond dimensions, the energy benefit of the Full Update scheme can become substantial. A lighter version of this approach is the implementation of a 3d extension of the neighborhood tensor update from Ref. [72], which combines higher accuracy than Simple Update with a still reasonable cost even in 3d. Another possible algorithmic development concerns the calculation of operator averages. In this study, we employed only the Simple Update environments to optimize bPEPS tensors; however, it is also possible to use Full Update. It may also be interesting to compare our observable calculations with the loop expansion modification of belief propagation from Ref. [73]. Furthermore, one can also investigate issues with potential CTMRG ill-convergence and compare the observable calculations from the rotated tensor networks to analyze algorithmic convergence.

The authors acknowledge support from the National Research Foundation of Ukraine under the call “Excellent science in Ukraine 2024–2026, project No. 0124U004372.

References↩︎