The study of strongly correlated electron systems has been a central theme in condensed matter physics for decades. These systems are characterized by interactions among electrons that cannot be treated in a perturbative manner, leading to emergent phenomena such as Mott insulating behavior [1], unconventional superconductivity [2, 3], quantum spin liquids [4, 5], and topological phases of matter [6, 7]. A particularly challenging frontier lies in understanding these phenomena in two-dimensional (2D) systems, where quantum fluctuations, geometric frustration, and topology interplay to create a rich tapestry of exotic phases.

Introduced by Steven R. White in 1992 [8], the Density Matrix Renormalization Group (DMRG) method has emerged as a cornerstone of numerical methods for studying strongly correlated systems. Originally designed for one-dimensional (1D) systems, DMRG achieves remarkable accuracy by optimizing wave functions that take the form of matrix product states (MPSs) [9,10,11,12,13]. Its success hinges on the area-law scaling of entanglement entropy between two subsystems [14, 15]: the entanglement entropy for gapped ground states grows proportionally to the boundary size, which remains constant in 1D systems. Thus, quantum states can be efficiently represented using MPS with a small bond dimension [14, 15], which has established DMRG as the benchmark method for solving 1D models, such the antiferromagnetic Heisenberg chains [16] and the Hubbard chain [17].

It is natural to apply DMRG to quasi-1D systems, such as ladder and cylinder systems [18]. However, DMRG exhibits sharply distinct performance in 1D and 2D systems, which originates from the different entanglement scaling with respect to spatial dimensionality. Unlike in 1D systems, the boundary size of the subsystem in 2D systems is not constant, and thereby the entanglement entropy grows with increasing system size. Consequently, the bond dimension of MPS required to faithfully capture a quantum state in 2D systems can grow exponentially with the system size, rendering the traditional DMRG method increasingly computational cost [12, 18]. Moreover, the convergence of DMRG is not guaranteed due to the presence of the notorious “local minima” in the energy landscape [19]. Despite these challenges, DMRG has been moderately successfully applied to narrow-width cylinders, providing important insights into 2D systems such as the kagome-lattice Heisenberg model [20,21,22] and square-lattice Hubbard model [23, 24].

To address the performance limitations of DMRG in two-dimensional systems, hybrid approaches have been developed to incorporate physically motivated variational Ansätze as initial inputs for DMRG simulations. Among these, the integration of Gutzwiller-projected wave functions with DMRG turns out to be a promising direction. The Gutzwiller projection, first introduced by Martin Gutzwiller [25], is a variational technique designed to eliminate unphysical configurations of a trial wave function, namely doubly occupied sites that are prevented by an infinitely large Hubbard U^{Footnote 1}. Gutzwiller-projected wave functions have traditionally been combined with variational Monte Carlo (VMC) methods [26, 27], serving as trial states to investigate quantum spin liquids and high-T\(_c\) superconductors [26]. On the other side, recent methodological progresses have shown that these trial wave functions can be formulated within the framework of tensor network states [28,29,30,31,32,33]. Specifically, converting them into MPS allows DMRG simulations to leverage high-quality initial states, guiding this optimization process toward a physically relevant state. This approach helps DMRG avoid local minima and accelerates convergence [19]. In addition to focusing on energetics, the quality of the trial state can be directly assessed by evaluating the wave function fidelity between itself and the DMRG-optimized state. This hybrid approach, known as “Gutzwiller-guided DMRG,” combines the computational power of DMRG with the physical insight encoded in the variational states, emerging as a promising tool for studying 2D systems.

Recent studies have demonstrated the potential of this hybrid method to address long-standing questions in 2D strongly correlated systems. For example, the spin-1/2 antiferromagnetic Heisenberg model on the kagome lattice, a paradigmatic system for quantum spin liquids, has been studied using this approach to reveal strong evidence for a chiral spin liquid (CSL) phase [34]. Remarkably, the boundary semion excitations in this CSL phase can be easily constructed without altering the physical Hamiltonian, highlighting the efficacy of this method in probing topological order. Similarly, the SU(4) Kugel-Khomskii model on the triangular lattice has been identified as hosting a nematic spin-orbital liquid state with an emergent parton Fermi surface [35], and the same model on the honeycomb lattice is found to be a \(\pi\)-flux state with Dirac-type excitations [36]. The Kitaev honeycomb model, known for its exact solution of the entire energy spectrum, has also been investigated using this method, providing valuable insights into its topological ground-state degeneracy [19].

In this review, we provide a comprehensive overview of the so-called Gutzwiller-guided DMRG method. We begin by discussing the theoretical framework of DMRG and Gutzwiller projection, followed by a detailed description of their integration. We then highlight applications of the method to prominent 2D models, exploring its role in uncovering exotic phases and topological orders. Finally, we discuss the challenges and future directions for this emerging approach, emphasizing its potential to revolutionize the study of 2D strongly correlated systems.

Owing to competing interactions and quantum fluctuations, strongly correlated electron systems are notoriously challenging to study using conventional techniques such as mean-field theory or perturbation theory. To effectively deal with these systems, theoretical and computational approaches usually rely on variational Ansätze to capture the essential entanglement features of the actual ground state for a given Hamiltonian. In this section, we outline the theoretical framework of Gutzwiller-projected wave functions and the DMRG method, followed by a discussion of their integration within a unified hybrid framework.

2.1 Gutzwiller-projected wave functions

The Gutzwiller projection, proposed by Martin Gutzwiller in the 1960s [25], provides a way to take into account correlations on top of non-interacting fermions. Broadly speaking, Gutzwiller-projected fermionic wave functions take the following form:

This section outlines the key computational steps of “Gutzwiller-guided DMRG,” including the construction of Gutzwiller wave functions, conversion into MPS, and the integration into the DMRG framework.

3.1 Construction of Gutzwiller wave functions

The process begins with the construction of a suitable Gutzwiller-projected wave function. The unprojected wave function, \(|\Psi _0\rangle\), is typically chosen as a mean-field solution to a quadratic Hamiltonian. Common choices for \(|\Psi _0\rangle\) include BCS-like states for paired fermions, Slater determinants for free fermions, or parton wave functions for spin systems [26, 28]. These states encode the non-interacting physics of the system and serve as a foundation for introducing strong correlations.

Let us elaborate on this aspect using the example of \(S=1/2\) spin systems with SU(2) symmetry. To construct a fermionic mean-field theory, we first express \(S=1/2\) spin operators at site i in terms of two-component fermionic partons \(\psi _i = (f_{i\uparrow }, f_{i\downarrow })^T\):

By leveraging physically motivated Gutzwiller-projected wave functions as initial ansatz, this hybrid method addresses key challenges in 2D strongly correlated systems, as evidenced by the improved computational efficiency and precision by orders of magnitude (see Table 1). In this section, we discuss its applications to several famous 2D systems, illustrating its effectiveness in uncovering ground-state properties, quasiparticle excitations, and topological order.

4.1 Kitaev honeycomb model

The \(S=1/2\) Kitaev honeycomb model [7] is a paradigmatic system in the study of quantum spin liquids, offering an exactly solvable framework for exploring fractionalization and topological ground-state degeneracy. The Hamiltonian of this model is given by

\(\begin{aligned} H_{\text {Kitaev}} = J_x \sum \limits _{\langle ij \rangle _x} \sigma _i^x \sigma _j^x +J_y \sum \limits _{\langle ij \rangle _y} \sigma _i^y \sigma _j^y +J_z \sum \limits _{\langle ij \rangle _z} \sigma _i^z \sigma _j^z, \end{aligned}\)

(12)

where \(\sigma _i^a\) (\(a = x, y, z\)) are three Pauli matrices, and \(\langle ij \rangle _a\) denotes nearest-neighbor bonds of type a (see Fig. 1). The bond-dependent anisotropic interactions lead to a highly frustrated system, where the ground state is known to be a Kitaev spin liquid with emergent Majorana fermions coupled to a static \(\mathbb {Z}_2\) gauge field.

Kitaev honeycomb model on a cylindrical geometry, in which the x-boundary is open while the y-boundary remains periodic. Black dots and white circles stand for A and B sublattices. The a-type (\(a=x,y,z\)) Ising interactions are defined on the a type of nearest-neighbor bonds. The purple zigzag line indicates a closed loop \(\mathcal {C}\) along which the Wilson loop operator Wy is defined (see Eq. (14)). Figure from Ref. [19]

 Full size image 

The exact solvability of the Kitaev honeycomb model provides an ideal benchmark for testing and validating the Gutzwiller-guided DMRG method. This observation is clearly illustrated by using Kitaev’s four-Majorana representation of spin as \(\sigma ^{a}_j=ic^{a}_jc^{0}_j\), where \(c^a\) and \(c^0\) are called gauge and itinerant Majorana fermions, respectively. Note that this representation is equivalent to the parton representation in Eq. (2), up to an SU(2) gauge transformation. Nevertheless, the local single-occupancy constraint becomes \(D_j\equiv c^{x}_jc^{y}_jc^{z}_jc^0_j=1\), and thereby the Gutzwiller projector for restoring the physical Hilbert space reads \(P_G = \prod \nolimits _j \frac{1+D_j}{2}\). Under this representation, \(H_{\text {Kitaev}}\) becomes

\(\begin{aligned} H_{\text {Kitaev}} = -i \sum \limits _a\sum \limits _{\langle jk \rangle \in a}J_a{u}_{jk} c^0_j{}c^0_k, \end{aligned}\)

(13)

where the static \(\mathbb {Z}_2\) gauge field \(u_{jk} \equiv ic_{j}^{a}c^{a}_{k}\) lives on an a-type bond. Noting that this gauge field commutes with the Hamiltonian for all different bonds, one can replace the operators of \(u_{jk}\) by their eigenvalues \(\pm 1\). After fixing the static \(\mathbb {Z}_2\) gauge field, \(H_{\text {Kitaev}}\) becomes a quadratic Hamiltonian of itinerant Majorana fermions, whose ground state is a BCS wave function.

4.1.1 Benchmark results

As a Gutzwiller-projected BCS wave function, the ground state of the Kitaev honeycomb model can be converted into an MPS using the MPO-MPS method [19]. On a finite cylinder with \(L_y\times {}L_x=4\times 10\), the MPO-MPS method leads to a tiny accumulated truncation error of \(\epsilon _{\text {trunc}}\sim {}10^{-2}\) with a small bond dimension \(\tilde{D}=200\).

We further perform DMRG optimization with the MPSs prepared from Gutzwiller-projected states, gradually increasing the bond dimension from \(\tilde{D}\) to D. For comparison, we also generate an MPS with random tensors with bond dimension \(\tilde{D}\), and then optimize it with the same DMRG scheme until a converged MPS of bond dimension D is obtained. The calculations are performed on an \(L_x \times L_y =6 \times 6\) cylinder. The truncation errors of DMRG are kept to smaller than \(10^{-9}\) during the optimization of DMRG. For comparison, we calculate the relative energy deviation \(\delta E_g\) between the DMRG energy and exact ground-state energy. We find that (1) the converged MPS reveals that Gutzwiller Ansätze achieve the optimal result with a final bond dimension of \(D=6500\), notably smaller than the required \(D=8000\) for random MPS initialization; (2) the required number of sweeps for Gutzwiller Ansätze is roughly half of those for random MPSs; and (3) the relative energy deviation \(\delta E_g\) for the Gutzwiller Ansätze is two orders of magnitude smaller than that of a random MPS. In general, this benchmark study fully demonstrates the superiority of Gutzwiller-guided DMRG.

4.1.2 Distinguishing topological sectors

The Kitaev honeycomb model exhibits topological degeneracy in manifolds with non-trivial topology, such as cylinders or tori. Here, we adopt cylindrical boundary conditions, where the honeycomb lattice is embedded on a finite cylinder with \(L_x\) (\(L_y\)) unit cells along the open (periodic) directions (see Fig. 1). The Hamiltonian \(H_{\text {Kitaev}}\) now commutes with the Wilson loop operator wrapping around the cylinder, e.g.,

\(\begin{aligned} W_y = -\prod \limits _{j \in \mathcal {C}} \sigma ^y_j \end{aligned}\)

(14)

with \(\mathcal {C}\) being a closed loop as shown in Fig. 1. The eigenvalue of \(W_y\) is just the product of the static \(\mathbb {Z}_2\) gauge fields along the loop,

\(\begin{aligned} W_y\vert \Psi _G\rangle = \Phi _y\vert \Psi _G\rangle \end{aligned}\)

(15)

with \(\Phi _y=\prod \nolimits _{\langle jk \rangle \in \mathcal {C}} u_{jk} = \pm 1\) labeling two different topological ground-state sectors.

The Gutzwiller-guided DMRG method is an ideal tool to target ground states in distinct topological sectors [19]. Moreover, during the DMRG sweeps initialized with Gutzwiller-projected states, the eigenvalue of the Wilson loop operator (\(\Phi _y=\pm 1\)) is preserved, i.e., the MPS stays in the respective sector. This is very useful for studying topologically ordered states with (quasi-)degenerate ground states on cylinders. In contrast, we found that the DMRG initialized with a random MPS always converges to an MPS in \(\Phi _y=-1\) sector on a \(6 \times 6\) cylinder. However, exact results indicate that for a finite cylinder, the ground-state energy in the \(\Phi _y=-1\) sector is higher than that in the \(\Phi _y=1\) sector. For instance, the energy difference on the \(6 \times 6\) cylinder is given by \(E_g(\Phi _y=-1)-E_g(\Phi _y=1)\approx 0.084\).

The above result implies that the DMRG with a random initial ansatz gets stuck in a local minimum. On contrast, when the initial state already captures key features of the expected ground state, e.g., topological properties or entanglement structures, the DMRG optimization is guided toward a more relevant region of the energy landscape, reducing the likelihood of converging to unphysical local minima. It is noted that any initial state introduces bias, including the random MPS, which could potentially limit the exploration of the full Hilbert space. However, as demonstrated in this example, in the context of DMRG for 2D systems, this controlled bias is advantageous, as it leverages prior knowledge to enhance efficiency and accuracy, much like informed initial guesses in other variational methods (e.g., in Monte Carlo simulations). This is because the Gutzwiller ansatz provides a “warm start” that aligns with the system’s topological and entanglement features, thereby narrowing the search space and mitigating the risk of local minima that arise from the complex, high-dimensional optimization in DMRG.

4.2 SU(4) Kugel-Khomskii model

The SU(4) Kugel-Khomskii model extends the physics of spin-orbital interactions beyond the conventional SU(2) framework by incorporating orbital degeneracies [59, 60]. In this model, each lattice site hosts both spin and orbital degrees of freedom, and their interactions are described by a Hamiltonian that respects an enlarged SU(4) symmetry:

\(\begin{aligned} H = \frac{1}{2}\sum \limits _{\langle i j \rangle } \left( \boldsymbol{\sigma }_i \cdot \boldsymbol{\sigma }_j +1 \right) \left( \boldsymbol{\tau }_i \cdot \boldsymbol{\tau }_j + 1 \right) , \end{aligned}\)

(16)

where \(\langle ij\rangle\) denotes a nearest-neighbor bond. Here, \(\boldsymbol{\sigma }\) and \(\boldsymbol{\tau }\) represent Pauli matrices for spin and orbital degrees of freedom, respectively. In particular, the equal footing spin and orbital degrees of freedom play a symmetric role as well and are promoted to SU(4) symmetry rather than the usual \(\mathrm {SU(2)} \times \mathrm {SU(2)} \simeq \mathrm {SO(4)}\) symmetry at this fine-tuned point.

Generally, the larger symmetry amplifies quantum fluctuations and might stabilize spin-orbital liquid ground states. Therefore, our Gutzwiller-guided DMRG method naturally provides a versatile tool for investigating possible emergent gauge fields and fractionalized excitations in the SU(4) Kugel-Khomskii model.

We first briefly review the fermionic parton construction for SU(4) quantum magnets. This starts by introducing four-flavor fermionic partons \(f^\dagger _{i,m} \; (m=1,\ldots ,4)\) and express the spin and orbital operators as

\(\begin{aligned} & \sigma ^a_{i} \rightarrow \boldsymbol{f}^\dagger _i\sigma ^a\tau ^0\boldsymbol{f}_i,\nonumber \\ & \tau ^b_{i}\rightarrow \boldsymbol{f}^\dagger _i\sigma ^0\tau ^b\boldsymbol{f}_i,\nonumber \\ & (\sigma ^a\tau ^b)_i \rightarrow \boldsymbol{f}^\dagger _i\sigma ^a\tau ^b\boldsymbol{f}_i, \end{aligned}\)

(17)

where \(\sigma ^0\) (\(\tau ^0\)) is the identity matrix in spin (orbital) space, and \(\boldsymbol{f}^\dagger _i=(f^\dagger _{i,1},\cdots ,f^\dagger _{i,4})\) are fermion creation operators. In order to restore the original Hilbert space, a single-occupancy constraint, \(\sum \nolimits _{m=1}^4 f^\dagger _{i,m}f_{i,m}=1\), has to be imposed at each site, which also defines the Gutzwiller projection for this SU(4) system.

With the help of the parton representation, a mean-field decomposition can be carried out to obtain various effective quadratic Hamiltonian of partons, which determines mean-field ground states and low-energy excitations. In the study of the SU(4) Kugel-Khomskii model on the triangular lattice, we consider a stripy state that breaks one of the lattice translational symmetries by doubling the unit cell, as shown in Fig. 2a. Without loss of generality, we take x-direction as the symmetry breaking direction, and write down the mean-field Hamiltonian for such a stripy state as [35]

\(\begin{aligned} H_{\text {sp}} =-\sum \limits _{\langle ij \rangle ,m} \left[ 1+(-1)^{r_i}\delta \right] \left( f^\dagger _{i,m}f_{j,m} + f^\dagger _{j,m}f_{i,m} \right) , \end{aligned}\)

(18)

where \(r_i\) is the column index of lattice site i, \(r_{j}\ge {}r_{i}\) is assumed for nearest-neighbor bond \(\langle {}ij\rangle\), and \(\delta {}\) denotes the stripy strength. A finite \(\delta\) deforms the circular Fermi surface at \(\delta =0\). We find that the deformed Fermi surface consists of open orbits in reciprocal space when \(\delta >\delta _{c}\approx 0.09\), as shown in Fig. 2b.

a A sketch for the stripy parton state on a \(L_y=4\) cylinder. The black and gray bonds represent \(1+\delta\) and \(1-\delta\) nearest-neighbor bonds in Eq. (18), respectively. b The single-particle spectrum for the lower band of the stripy parton state with \(\delta =0.15\) in the (unfolded) first Brillouin zone. The blue curve stands for the deformed Fermi surface at 1/4 filling, which consists of open orbits in the reciprocal space. Dashed lines enclose the folded Brillouin zone. Green lines represent momenta on XC4 cylinders that are allowed by the APBC along y direction for partons. Figure from Ref. [19]

 Full size image 

This parton mean-field Hamiltonian (18) gives rise to a Fermi sea ground state at quarter filling, whose Gutzwiller projected version can be converted into MPS with the MPO-MPS method [28]. To perform MPS-related calculations on 2D lattices, we must must place the system on cylindrical geometries and work with finite-width cylinders. Note that fermionic partons are coupled to an emergent U(1) gauge field which can lead to a global gauge flux, \(\Phi\), through the cylinder along the periodic boundary conditions. The time-reversal symmetry only allows \(\Phi =0\) or \(\Phi =\pi\), corresponding to periodic boundary condition (PBC) or anti-periodic boundary condition (APBC) for fermionic partons, respectively. Note that \(\Phi\) is not a physical observable and is well-defined only with respect to a parton Hamiltonian with fixed U(1) gauge. Therefore, choosing APBC along the y direction for fermionic partons still yields PBC for spin and orbital degrees of freedom.

For finite cylinders, choosing PBC and APBC for fermionic partons corresponds to different ways of cutting the Brillouin zone. Crucially, for a finite-width cylinder, the energy can be reduced by avoiding the cutting of momenta at any existing gapless points. Although it is difficult to prove this argument analytically, numerical validation is feasible using the Gutzwiller-guided DMRG method.

Using the combinative method, we obtain the ground states of the SU(4) Kugel-Khomskii model on a triangular lattice with different cylindrical boundary conditions and various system sizes. We then evaluate the variational ground-state energy of the stripy ansatz and find that it reaches its minimum at \(\delta \approx 0.15\), indicating a stripy spin liquid state. Moreover, we calculate the wave function fidelity defined by

\(\begin{aligned} F = \langle \Psi _{\text {Gutzwiller}}|\Psi _{\text {DMRG}}\rangle . \end{aligned}\)

(19)

Notably, F is also maximized at \(\delta {}\approx 0.15\) with a value of \(F\approx 0.85\) on a \(4\times 8\) cylinder, in good agreement with the energy comparison.

It is remarkable that the optimal value of \(\delta\) does not open an energy gap on the parton Fermi surface. Thereby, a deformed parton Fermi surface consisting of open orbits is expected at quarter-filling in the two-dimensional limit. However, due to the global emergent U(1)-flux \(\Phi\), none of the momenta allowed in these finite \(L_y\) cylinders cuts the Fermi surface. For instance, \(\Phi =\pi\) has a lower energy for \(L_y=4\) cylinders (which leads to a gapped parton Fermi surface), while \(\Phi =0\) is energetically favored for \(L_y=6\) cylinders. Similar results have been found in the SU(4) Kugel-Khomskii model on the honeycomb lattice, which supports a Dirac spin liquid [36, 61]. On finite cylinders, the Dirac cones are also gapped by the emergent global U(1)-flux \(\Phi\) under the same mechanism.

The above results indicate another important advantage of the Gutzwiller-guided DMRG method, namely that it provides a diagnosis tool for analyzing the quality of Gutzwiller ansatz for a given Hamiltonian. As discussed above, the SU(4) Kugel-Khomskii model on the triangular and honeycomb lattices, while fundamentally gapless, can exhibit gapped features on finite cylindrical geometries, which might lead to misidentification as gapped states (if no special care has been taken). Consequently, one may misidentify the SU(4) Kugel-Khomskii model on quasi-1D geometries as a gapped QSL or other states when focusing only on local quantities like energetics and (short-distance) correlation functions. On the other hand, the Gutzwiller-guided DMRG method serves as a powerful tool for evaluating wave function fidelities F. Unlike local correlation functions, a high value of fidelity F, encoding non-local information, provides a more reliable indicator, effectively distinguishing true bulk behavior from finite-size artifacts.

4.3 Kagome Heisenberg antiferromagnet

The \(S=1/2\) kagome Heisenberg antiferromagnet (KHAF) is one of the most intriguing and challenging systems in the field of quantum magnetism, renowned for its potential to host exotic quantum spin liquids. The model describes a system of spin-\(\frac{1}{2}\) particles arranged on a two-dimensional kagome lattice, where the spins interact via antiferromagnetic nearest-neighbor couplings. The Hamiltonian is given by

\(\begin{aligned} H_{\text {KHAF}} = \sum \limits _{\langle ij \rangle _1} \textbf{S}_i \cdot \textbf{S}_j, \end{aligned}\)

(20)

where \(\langle ij \rangle _1\) denotes the 1st nearest-neighbor bonds. The geometric frustration inherent in the kagome lattice, combined with quantum fluctuations, suppresses conventional magnetic ordering. This makes the KHAF a prime candidate for realizing quantum spin liquid states.

Over decades, the ground state of the KHAF has been studied with many different approaches, but its nature is still under debate. Early theoretical studies suggested \(\sqrt{3}\times {}\sqrt{3}\) magnetic order [62] or valence-bond crystal with spontaneously broken lattice symmetry [63,64,65,66,67]. However, later numerical studies (in particular, DMRG and tensor networks) provided evidence that KHAF may instead host a quantum spin liquid ground state [20,21,22, 68,69,70,71,72]. The distinction between gapped and gapless spin liquid phases has been a focal point of these studies.

With the Gutzwiller-guided DMRG method at hand, we have revisited the KHAF in Ref. [34]. It provides strong numerical evidence that the ground state of the KHAF might be a chiral spin liquid (CSL), with spontaneously broken time-reversal symmetry (TRS). Previous works report that this Kalmeyer-Laughlin-type CSL is stabilized in an extended kagome Heisenberg model with the following Hamiltonian [73, 74]:

\(\begin{aligned} H_{\text {ex}} = H_{\text {KHAF}} + J'\left( \sum \limits _{\langle ij \rangle _2} \textbf{S}_i \cdot \textbf{S}_j+\sum \limits _{\langle ij \rangle _3} \textbf{S}_i \cdot \textbf{S}_j\right) , \end{aligned}\)

(21)

where \(\langle ij \rangle _n\) denotes the nth nearest-neighbor bonds (the third-neighbor couplings only exist within the hexagons). However, our results in Ref. [34] demonstrate that this CSL phase includes the \(J'=0\) point, i.e., the KHAF with only 1st nearest-neighbor interactions.

4.3.1 Chiral spin liquid phase and evidence

A CSL is a gapped quantum spin liquid characterized by broken time-reversal as well as non-trivial topological order. In the parton representation, a Kalmeyer-Laughlin type CSL can be constructed using the following fermionic parton Hamiltonian [27]:

\(\begin{aligned} H_{\text {CSL}} = \sum \limits _{\langle {}ij\rangle _1}\sum \limits _{\sigma =\uparrow ,\downarrow }(\chi _{ij}f^{\dagger }_{i\sigma }f_{j\sigma } + \text {H.c.}), \end{aligned}\)

(22)

where \(\chi _{ij}=e^{\phi _{ij}}\) are link variables on the 1st nearest-neighbor bonds. Due to SU(2) gauge redundancy [6], the ground states of Eq. (22) are distinguished by the gauge fluxes threading through elementary triangles and hexagons on the kagome lattice. The CSL state of particular interest corresponds to the ansatz with \(\frac{\pi }{2}\)-flux through triangles and zero-flux through hexagons.

Being a topologically ordered state, the degeneracy of the Gutzwiller-projected ground state, governed by \(H_{\text {CSL}}\) in Eq. (22), depends on the topology of the manifold. When placed on a cylinder, it supports four exact zero modes at the boundary, denoted as \(d_{L\sigma }^{\dagger }\) and \(d_{R\sigma }^{\dagger }\), where L and R represent the left and right boundary of the cylinder, respectively. Then, the minimally entangled states (anyon eigen basis) [75] of such a CSL are constructed as follows [28, 76]:

\(\begin{aligned} |\Psi _{1}\rangle = P_{G} d_{L\uparrow }^{\dagger }d_{L\downarrow }^{\dagger }|\Phi \rangle , \quad |\Psi _{2}\rangle = P_{G} d_{L\uparrow }^{\dagger }d_{R\downarrow }^{\dagger } |\Phi \rangle ~, \end{aligned}\)

(23)

where \(|\Phi \rangle\) is a Fermi sea state with all negative energy modes of Eq. (22) being occupied, and \(P_G\) is the Gutzwiller projector ensuring single occupancy. Using the MPO-MPS method [28], the CSL anyon eigenbasis in Eq. (23) can be efficiently converted into MPSs to carry out the Gutzwiller-guided DMRG process.

In the KHAF, we have identified the CSL phase by several key signatures:

1. Spontaneous TRS breaking at \(J'=0\): We perform DMRG calculations initialized with randomly generated MPSs (denoted as “Random-DMRG” below) as well as the Gutzwiller-guided DMRG initialized with \(|\Psi _1\rangle\). The ground-state energies obtained by Random-DMRG and Gutzwiller-guided DMRG are almost identical, and the ground-state energy per site is consistent with the best available DMRG results [20, 21]. With different initializations, the ground state \(|\Psi _{\text {I}}\rangle\) and its time-reversal partner \(|{\Psi _{\text {I}}}^{*}\rangle\) are obtained at \(J'=0\), indicating the breaking of TRS.

   The spin chirality operator, which is defined on elementary triangles \(\triangle\) of sites (i, j, k) as

   \(\begin{aligned} \chi _{ijk} = \textbf{S}_i \cdot (\textbf{S}_j \times \textbf{S}_k), \end{aligned}\)

   (24)

   
   

   serves as a marker of TRS breaking [77]. However, directly evaluating the spin chirality order is challenging in previous DMRG studies that predominantly employed real-number MPS. Because the spin chirality operator transforms as \(\chi _{ijk} \rightarrow -\chi _{ijk}\) under time reversal, the expectation value of \(\chi _{ijk}\) with respect to a real-number MPS is always zero. Instead, these previous works calculated the chiral-chiral correlation function, resulting in an expectation value of spin chirality that is too small to be distinguished from (numerical) zero. On the other hand, our Gutzwiller-guide DMRG can target the two-fold Kramers degeneracy, which allowed us to directly evaluate the local spin chirality \(\langle \chi _{ijk}\rangle\). We find that \(\langle \chi _{ijk}\rangle \sim {}10^{-4}\) is small but does not vanish numerically (note that the truncation error in our DMRG calculations is \(10^{-6}\)), and is stable against finite-size scaling.

2. Topological degeneracy: Another essential signal of the CSL at \(J' = 0\) is the existence of (quasi-)degenerate ground states, \(|\Psi _{S}\rangle\), in the semion sector. For the CSL state at \(J'=0.2\), \(|\Psi _{S}\rangle\) can be clearly established, as the characteristic level counting of the entanglement spectra is consistent with the chiral SU(2)\(_1\) Wess-Zumino-Witten conformal field theory. We further found that \(|\Psi _{S}\rangle\) can be adiabatically evolved from \(J' = 0.2\) to \(J' = 0\) using the Gutzwiller-guided DMRG method, in which the wave function fidelity F between two neighboring-\(J'\) states, e.g., \(|\Psi _{S}(J')\rangle\) and \(|\Psi _{S}(J'-0.01)\rangle\), is always close to unity for \(0.0\le J' \le 0.2\). The entanglement spectra of \(|\Psi _{S}\rangle\) also show qualitative consistency in the region of \(0 \le J' \le 0.2\). Moreover, we showed that the per-site energy variances of \(|\Psi _{S}\rangle\) for the KHAF model at \(J'=0\), \(\sigma ^{2}/N\) is as small as \(\sim 10^{-4}\), suggesting that \(|\Psi _{S}\rangle\) is indeed an eigenstate of the Hamiltonian in Eq. (20).

   The identification of the CSL phase in the KHAF is robust against finite-size scaling. For instance, by pushing the bond dimension up to \(D=18000\) and the circumference of cylinder up to \(L_y=8\), we confirmed that the results on spin chirality and semion sectors are more evident as D and \(L_y\) increase. Nevertheless, challenges remain in fully characterizing the CSL phase in the thermodynamic limit. The finite circumference of the cylinder geometry limits the resolution of long-range correlations, and the growth of the computational cost with bond dimension imposes practical constraints. Future studies employing wider cylinders or alternative tensor network approaches, such as projected entangled pair states, may help address these limitations.

In this review, we summarize recent progress on integration of Gutzwiller-projected wave functions with the DMRG method. By combining the strengths of Gutzwiller projected fermionic parton wave functions and MPS techniques, we have devised a promising hybrid approach for investigating 2D strongly correlated systems. The key advantages are (1) the notorious “local minimum” issue in DMRG can be circumvented, and the precision of DMRG can be improved by orders of magnitude without extra computational cost; (2) apart from local correlation functions, the new method can be used to compute more essential non-local observables such as wave function fidelity, which allows us to directly diagnose the quality of the parton wave functions; and (3) it provides a versatile tool to target (topologically) degenerate ground states in different (topological) sectors.

A significant result of our investigation has been the accurate capture of complex phenomena such as topological order, fractionalized excitations, and the subtle interplay between competing interactions. For the cases of the Kitaev honeycomb model, the SU(4) Kugel-Khomskii model, and the KHAF model, the hybrid approach allows us to provide deeper insights into emergent behaviors in these systems. The flexibility of the method in adapting to various Hamiltonians and lattice geometries demonstrates its potential as a universal tool in the study of quantum spin liquids and other strongly correlated phases of matter.

The combination of Gutzwiller-projected wave functions with DMRG represents a significant step forward in the numerical study of many-body physics, providing a pathway to systematically explore and characterize quantum phases that emerge from strong electron correlations. This hybrid approach is expected to play a crucial role in resolving future challenges in the field of quantum magnetism and beyond.

Along this direction, there are several key issues that warrant further investigation:

1. Developing more efficient algorithm for converting fermionic Gaussian states into MPSs. The existing method introduced in Ref. [28] is effective yet somewhat brute-force and might not constitute the best strategy, especially within 2D systems. A potential alternative approach is to divide the entire system into several separate subsystems as small as possible. In this context, the Wannier orbitals correspond precisely to the Schmidt vectors of the reduced density matrices for each subsystem, which can be easily computed using the fermionic Gaussian state theory.

2. The transformation of the Gutzwiller-projected parton state into a 2D tensor network state, such as projected entangled pair state (PEPS), holds great promise. Compared to MPS, the optimization of PEPS relies even more heavily on the choice of initial states. Considering PEPS can capture entanglement entropy beyond the area law, this transformation could potentially unlock new ways of understanding and simulating complex quantum systems [32, 33].

3. The Gutzwiller projection is a powerful tool, but its application has been mainly restricted to the case of an infinitely large Hubbard \(U\rightarrow \infty\). However, the effects of charge fluctuations are inherently unavoidable in realistic systems. Extending this framework to accommodate finite Hubbard U interactions [78] and/or deviations from half-filling is a pressing issue for future research.

4. The Gutzwiller-guided DMRG has been demonstrated to be highly efficient in targeting the ground state in specific topological sectors. Building on this, the next step involves identifying potential applications for a range of topological systems, including but not limited to symmetry-protected topological states [79] and symmetry-enriched topological states.