Quantum spin liquid with emergent chiral order in the triangular-lattice Hubbard model

The interplay between spin frustration and charge fluctuation gives rise to an exotic quantum state in the intermediate-interaction regime of the half-filled triangular-lattice Hubbard model, while the nature of the state is under debate. Using the density matrix renormalization group with $\mathrm{SU}{\left(2\right)}_{\mathrm{spin}}\otimes \mathrm{U}{\left(1\right)}_{\mathrm{charge}}$ symmetries implemented, we study the triangular-lattice Hubbard model defined on the long cylinder geometry up to circumference $W=6$. A gapped quantum spin liquid, with on-site interaction $9\lesssim U/t\lesssim 10.75$, is identified between the metallic and the antiferromagnetic Mott insulating phases. In particular, we find that this spin liquid develops a robust long-range spin scalar-chiral correlation as the system length $L$ increases, which unambiguously unveils the spontaneous time-reversal symmetry breaking. In addition, the degeneracy of the entanglement spectrum supports symmetry fractionalization and spinon edge modes in the obtained ground state. The possible origin of chiral order in this intermediate spin liquid and its relation to the rotonlike excitations have also been discussed.

Figure 1

Model and phase diagram. (a) Triangular-lattice cylinder with open (periodic) boundary conditions along the $x$ ($y$) direction and lattice spacing $a=1$. Having a straight open edge, i.e., YC4 geometry, the columnar spacing is $a^{\prime }=\sqrt{3}/2$. (For later reference, an XC geometry has straight horizontal chains at distance $a^{\prime }$ with a zigzag open boundary left and right.) The DMRG simulations employ a 1D “zigzag” mapping as indicated by the gray-shaded path (with part of the site ordering also indicated). The yellow star denotes the central reference site for which the spin-spin correlation is calculated, while the pink triangles for evaluating the chiral correlations are ordered symmetrically away from the center. The arrows on the bonds denote the current directions of the chiral order, having three colors for the three different directions. (b) The phase diagram of the TLU model consists of a metallic phase, a fully gapped CSL phase, and a ${120}^{\circ }$ spin-ordered phase, with the three colors of the sites denoting the three-sublattice structure. (c)–(e) Typical static spin structure factors $S_q$ for $U=8,10$, and 12 in the three phases.

Figure 2

Hubbard $U$ dependence of different quantities on the $\mathrm{YC}4\times L$ cylinders. (a) Spin structure factors $S_q$ at the $K$ and $M$ points. The $L=64$ data are also shown at $U=10$. The asterisks show data points for the significantly longer system that exhibits long-range chiral correlations as analyzed in Fig. 4. (b) Double occupancy shown as $n_d{U}^2$. This also includes XC4 data, which exhibit only one discontinuity around the metal-insulator transition. (c) Distribution of $n_q$ vs $q_x$ in the momentum space along $q_y=\pi /2$. The dashed line plots $n_q$ for $U=0$ on the torus system with the same size $4\times 18$. The white center region marks the first Brillouin zone.

Figure 3

Single-particle Green's function, spin correlation, and entanglement entropy obtained on the $\mathrm{YC}4\times 64$ cylinder at $U=10$. Log-linear plot of (a) single-particle Green's function $|\langle c_i^{†}{c}_j^{}\rangle |$, and (b) spin correlation $|\langle {\mathit{S}}_i·{\mathit{S}}_j\rangle |$ as a function of distance $d_{ij}$, which show exponential decay with short decay lengths ${\xi }_c\simeq 0.81$ and ${\xi }_s\simeq 1.66$, respectively. (c) Bipartite entanglement entropy $S_E=-{\sum }_i{\rho }_{i}ln{\rho }_i$ for the reduced density matrix $\rho$ when cutting the system at bond $l_x$ vs. block size $l_x$. This shows a well-developed plateau for $12\lesssim l_x\lesssim 52$, and hence obeys the area law in the bulk. We include a linear extrapolation $1/D^{*}\to 0$ (black symbols). $S_E$ is also obtained from a complex wave function (red symbols; all other data for a real wave function) with the value of the plateau in the center reduced by $ln(2\pm 0.1)$.

Figure 4

Chiral correlation and entanglement spectrum of the chiral spin liquid state. (a) Chiral correlations $\langle {\chi }_i{\chi }_j\rangle$ for $U=10$ on the $\mathrm{YC}4$ cylinders with different system lengths $L=16,18,64$. The orange horizontal line indicates the value of $\langle {\chi }_i{\chi }_j\rangle \simeq 0.128$ for $d_{ij}=30$. (b) Chiral correlations for $U=9,10,11$ on the $\mathrm{YC}4\times 64$ cylinder obtained by keeping the bond dimensions up to $D^{*}=4096$ multiplets. Entanglement spectra of the $\mathrm{YC}4\times 64$ for both (c) real and (d) complex wave functions, grouped by charge sectors (for even $Q$), with spin labels color-coded as specified in the legend. The bars and respective numbers with the $Q=0$ column in (c) indicate group degeneracy. The subtracted ground levels for the two cases are ${\lambda }_0\simeq 0.08$ and ${\lambda }_0\simeq 0.15$, where a relative factor of $\sim 2$ is observed.

Figure 5

DMRG simulations (real-valued) on $\mathrm{YC}6\times 64$ at $U=9$. (a) Convergence of ground-state energy vs discarded weight $\delta \rho$. The lowest energies are extrapolated in a linear fashion towards $\delta \rho \to 0$, with the resulting ground-state energy per site as shown. (b) Maximum entanglement entropy around the system center vs $1/D^{*}$, with $D^{*}$ the number of multiplets kept in the simulation. The inset shows the discarded weight vs $1/D^{*}$. The lines in panels (a) and (b) combine data from equivalent sweeps, such as the two-site update (line with symbol) when increasing $D^{*}$, or the subsequent four bond updates [for all other lines, see the legend to (a)]. (c) Chiral correlations vs distance. While these data show exponential decay, the correlation length is strongly dependent on $D^{*}$ still. Here the various lines are derived from different combinations of up and down triangles vs distance. Different levels of color intensity refer to different $D^{*}$ ($D$) as indicated in the legend. (d) Analysis of the parameters from the exponential fits in (c) vs discarded weight $\delta \rho$.

Figure 6

DMRG simulation (complex-valued) on $\mathrm{YC}6\times 64$ at $U=9$. (a)–(c) Identical analysis as in Fig. 5. (d) Having complex arithmetic, this permits a nonzero expectation value for the chiral order parameter $\langle {\chi }_i\rangle$, showing the square for direct comparison with (c). We note that $\langle {\chi }_i\rangle$ has the same sign for all triangles and, consistent with (c), converges towards $|\langle \chi \rangle |\approx 0.25$ (horizontal red dashed line). The data in (d) correspond to the last three entries in the legend in (c).

Figure 7

Correlations in the $\mathrm{YC}6\times 64$ system after the last sweep in the complex DMRG simulation at $D^{*}=8192$ in Fig. 6: (a) Spin $\langle {\mathbf{S}}_i·{\mathbf{S}}_j\rangle$ and charge correlations $\langle {\mathbf{c}}_i^{†}·{\mathbf{c}}_j\rangle$ (where the dot product sums over spin) relative to the system center. The exponential fit $\propto e^{-x/\xi }$ (straight line in matched light color) yields the corresponding correlation length as shown with the legend. The correlations are computed along a tilted straight path of length $L-1$ in units of lattice spacing from left to right YC boundary [cf. Fig. 1]. (b) Chiral correlations, (c) spin correlations, (d) charge correlations, all at constant vertical distance $d_y=3$ ($j=2,5$) vs position $x$ along the cylinder. Here $x$ is the distance from the boundary of the cylinder converted to column index $x/a^{\prime }\in [1,L]$ [cf. Fig. 1].

Figure 8

Symmetric tensor representation of (a) SU(2) spinor $S^q$ [Eq. (A1)], (b) $P$-operator, (c) the 1j-symbol, and (d) the scalar chirality operator order $\chi$ [cf. Eq. (A2b)].

Figure 9

Entanglement entropy $S_E$ (left panels) and ground-state energy (right panels) vs $1/D^{*}$ for YC4 for (a),(b) $U=8$; (c),(d) $U=10$; (e),(f) $U=12$ at $L=18$; and (g),(h) $U=10$ at $L=64$. Entanglement entropy $S_E$ labeled “Max” are the maximal entanglement among all bonds throughout the system, while the “Mid” ones are values measured at the center of the systems. The linear extrapolations with $1/D^{*}$ are shown in the inset of all right panels (b),(d),(f),(h). Although the calculations are very challenging in the metallic phase with $U=8$, the data have reached very good convergence for $U=10$ and 12.

Figure 10

Entanglement spectrum (ES) along the system and in between full columns of a DMRG scan along a $\mathrm{YC}4\times 100$ cylinder. For this we tune $U\in [8,12]$ linearly along the cylinder. We keep the bond dimensions up to $D^{*}=5793$ ($D\lesssim 17053$). The lines are color-coded according their symmetry sector $(S,Q)$ as indicated.

Figure 11

Electron density in the momentum space $n_q$ on the $\mathrm{YC}4\times 18$ systems with (a) $U=0$, (b) $U=8$, (c) $U=9$, and (d) $U=10$. Gray hexagons represent the boundary of first Brillouin zone, and the horizontal dotted lines represent the allowed momenta $q_y$ for the YC4 cylinder.

Figure 12

Charge gap ${\mathrm{\Delta }}_{\mathrm{C}}$ calculated on the YC4 cylinders with different lengths $L=8,12,16,18$. The linearly extrapolated value with $L\to \infty$ is ${\mathrm{\Delta }}_{\mathrm{C}}\simeq 0.88$.

Figure 13

Energy difference ${\mathrm{\Delta }}_S$ between the total spin-0 and total spin-1 sectors. We calculate the YC4 cylinders with system length $L=8,12,16,18$. This shows an exponential decay with $\xi \simeq 5.8$ with increasing $L$ as indicated.

Figure 14

Entanglement spectrum calculated at the center of YC4 systems with (a) $L=18$ and (b) $L=64$ [cf. Fig. 4 in the main text], grouped by charge sectors (showing even $Q$ only), with spin labels color-coded as specified in the legend.

Figure 15

Chiral correlations calculated in a $U=9.5 \mathrm{XC}4\times 64$ system with bond dimensions $D^{*}=4096,5793,8192$ show no long-range correlations. The black line depicts the results of linear extrapolation $1/D^{*}\to 0$ from the above bond dimensions. Here, the filled symbol indicated the positive sign of chiral correlation, and otherwise the sign is negative.