Spin-Orbital Density Wave and a Mott Insulator in a Two-Orbital Hubbard Model on a Honeycomb Lattice

Inspired by the recent discovery of correlated insulating states in twisted bilayer graphene, we study a two-orbital Hubbard model on the honeycomb lattice with two electrons per unit cell. Based on the real-space density matrix renormalization group simulation, we identify a metal-insulator transition around $U_c/t=2.5–3$. In the vicinity of $U_c$, we find strong spin-orbital density wave fluctuations at commensurate wave vectors, accompanied by weaker incommensurate charge density wave fluctuations. The spin-orbital density wave fluctuations are enhanced with increasing system sizes, suggesting the possible emergence of long-range order in the two-dimensional limit. At larger $U$, our calculations indicate a possible nonmagnetic Mott insulator phase without spin or orbital polarization. Our findings offer new insight into correlated electron phenomena in twisted bilayer graphene and other multiorbital honeycomb materials.

Figure 1

The average of the double occupancy $n_d$ as a function of $U/t$ for torus (a) and cylinder (b) geometries. The insets show the $n_{d}U$ as a function of $U/t$. The shadow bar indicates the MIT.

Figure 2

The momentum distribution $n(\mathbf{k})$ for the metallic phase at $U/t=1$ (a) and insulating phase at $U/t=4$ (b). The black dots represent the momentum points we can access in the Brillouin zone of the finite size lattice; the contour plot is created by using triangulation interpolation. The cut of $n(\mathbf{k})$ along the line crossing $\mathrm{\Gamma }$ point for different ratios of $U/t$ are shown in (c), where the slope of $n(\mathbf{k})$ decreases with the increase of $U/t$ near the Fermi surface. Here, we consider the cylinders with $N_0=12\times 3\times 2$. (d) The same cut of $n(\mathbf{k})$ as (c) at $U/t=1$ and $U/t=4$ for different system sizes.

Figure 3

The contour plot of the structure factor $N_{\mathbf{q}}$ for $L_y=3$ cylinders at $U/t=2$ (a), $U/t=2.5$ (b), $U/t=3$ (c), and $U/t=4$ (d). $N_{\mathbf{q}}$ are featureless for both metallic phase (a) and insulating phase (d), while the CDW fluctuations arise in the vicinity of the MIT (c). The splitting of the CDW peaks around $M$ point signifies the incommensurate CDW instability. Here, we consider the cylinders with $L_x=12$ unit cells; the black dots represent the momentum points we can access in the Brillouin zone of the finite size lattice. The contour plot is created by using triangulation interpolation.

Figure 4

The contour plot of the static spin structure factor $S_{\mathbf{q}}$ for $L_y=3$ cylinders at $U/t=2$ (a), $U/t=2.5$ (b), $U/t=3$ (c), and $U/t=4$ (d). $S_{\mathbf{q}}$ are featureless for both metallic phase (a) and insulating phase (d), while the spin-orbital density wave fluctuations arise in the vicinity of the MIT point (b),(c). The peaks of $S_{\mathbf{q}}$ locate at $K$ or $M$ points indicating the commensurate spin-orbital density wave instability. Here, we consider the cylinders with $L_x=12$ unit cells; the black dots represent the momentum points we can access in the Brillouin zone of the finite size lattice. The contour plot is created by using triangulation interpolation. The peaks in (c) is higher than (b), indicating the stronger fluctuations with wave vector at $M$ points than at $K$ points.

Figure 5

The structure factor $N_{\mathbf{q}}$ (a) and $S_{\mathbf{q}}$ (b) for $L_y=4$ systems at $U/t=3$. Here, we consider the cylinders with $L_x=12$, the black dots represent the momentum points we can access in the Brillouin zone of finite size lattice system. The contour plot is created by using triangulation interpolation.

Figure 6

The ground state energy $E_0$ per site as a function of $U/t$ for both unpolarized $P=0$ and polarized $P=1$ sectors. Here, we consider the systems with $L_x=12$ and different $L_y$.