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

TL;DR

This study uses DMRG simulations to explore a two-orbital Hubbard model on a honeycomb lattice, revealing a metal-insulator transition, density wave fluctuations, and a potential nonmagnetic Mott insulator phase relevant to twisted bilayer graphene.

Contribution

It provides the first detailed DMRG analysis of a two-orbital Hubbard model on a honeycomb lattice, uncovering novel density wave phenomena and a possible Mott insulator phase.

Findings

Metal-insulator transition at Uc/t≈2.5-3

Strong spin/orbital density wave fluctuations near transition

Possible nonmagnetic Mott insulator at larger U

Abstract

Inspired by recent discovery of correlated insulating states in twisted bilayer graphene (TBG), 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 (DMRG) simulation, we identify a metal-insulator transition around $U_c/t=2.5\sim3$. In the vicinity of $U_c$, we find strong spin/orbital density wave fluctuations at commensurate wavevectors, accompanied by weaker incommensurate charge density wave (CDW) 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 insights into correlated electron phenomena in twisted bilayer graphene and other multi-orbital honeycomb materials.

Figures

Equations

$$H$$

$$H_0$$

$$H_{\mathrm{int}}$$

Full text

Introduction.—The metal-insulator transition (MIT) driven by electron repulsion is one of the long-standing issues in condensed matter physics Mott ; Mott1990 ; Imada1998 . In particular, the systems with multiple orbits may harbor different exotic phases and novel phase transitions due to the interplay between orbit, spin and charge degrees of freedom.

Recently, correlated insulators were discovered at low temperature in twisted bilayer graphene (TBG) with an integer number of electrons per superlattice unit cell  Cao1 ; Cao2 ; Yankowitz2018 . These insulating states are found near the “magic” twist angle  Cao1 ; Cao2 or under pressure  Yankowitz2018 , where the lowest mini-band of TBG has a very narrow bandwidth. As a result of the reduced kinetic energy, the correlation effect becomes significant and most likely is the driving force for the MIT. The nature of the insulating states is now under intensive study  Fu2018a ; Yang ; Senthil ; Kivelson ; Vafek2 ; Pizarro ; Zhang ; PALee ; Xu ; Irkhin ; Rademaker ; Kuroki ; Fernandes ; Fu2018c ; You . A paradigmatic theoretical model for MIT is the Hubbard model and its various generalizations to include multiple orbits and extended interactions. For TBG, a two-orbital extended Hubbard model on the honeycomb lattice  Fu2018a ; Fu2018b was recently derived from symmetry-adapted maximally-localized Wannier states of the low-energy bands  Fu2018b ; Vafek . Since these low-energy bands are well separated from higher bands by a sufficiently large energy gap  Cao2016 ; Koshino2017 , this two-orbital extended Hubbard model is a reasonable and adequate starting point for studying correlated electron phenomena in TBG Fu2018b ; YuanFu . In this model, the two sublattices of the emergent honeycomb lattice correspond to Wannier states centered at AB and BA stacking regions, and the two orbits are associated with states from opposite valleys of graphene.

In a completely different context, the two-orbital $SU(4)$-symmetric Hubbard model on the honeycomb lattice was recently proposed for $d^{1}$ transition metal compounds with edge-sharing anion octahedra such as $\alpha$-ZrCl3 Yamada , where spin and orbital combined together describe the $j=3/2$ quartets that emerge in the strong spin-orbit-coupling limit.

In addition to its potential relevance to real materials, it is also of fundamental importance to study the physics beyond $SU(2)$. In contrast to large-$S$ representations of $SU(2)$, the quantum fluctuations in $SU(N)$ increase with $N$ and thus opening new realm to explore more intriguing physics. In this work, we study the two-orbital Hubbard model with $SU(4)$ symmetry on the honeycomb lattice at quarter filling, i.e., with two fermions per unit cell. Based on large-scale density matrix renormalization group White1992 (DMRG) simulations, we examine the nature of the ground state with the ratio of on-site Coulomb repulsion $U$ and bandwidth $t$. We identify a MIT occurs around $U_{c}/t=2.5\sim 3$ and a nonmagnetic Mott insulator in the large $U$ regime. Interestingly, in the vicinity of $U_{c}$, we find a small region where both spin/orbital density wave (SODW) and charge density wave (CDW) fluctuations become strong. The SODW fluctuations have commensurate wavevectors near $U_{c}$, and become enhanced when increasing system size, indicating long-range SODW order in the two dimensional (2D) limit. In contrast, the CDW fluctuations have incommensurate wavevectors and are unstable against increasing system size. The SODW is robust against perturbations including nearest neighbor interactions or lightly doping. We do not find any sign of spin/orbital polarization in the ground state.

Our findings of the nonmagnetic Mott insulator at large $U$ is consistent with previous studies of the $SU(4)$ Heisenberg model on honeycomb lattice Mila ; Penc , which is the effective Hamiltonian of our Hubbard model in $U/t\rightarrow\infty$ limit. Meanwhile, our finding of the SODW near MIT at intermediate $U$ is a new result. Our work shows that the phase diagram of $SU(4)$ Hubbard model is distinctly different from SU(2) case Meng2010 ; Chen2014 ; Sorella2012 ; Assaad2013 ; Hassan2013 ; Toldin2015 ; Otsuka2015 ; Shirakawa2015 ; Szasz2018 , where a magnetic ordered phase often appears in large $U$ limit, a direct phase transitionSorella2012 ; Assaad2013 ; Hassan2013 ; Toldin2015 ; Otsuka2015 or an intermediate phaseMeng2010 ; Chen2014 ; Shirakawa2015 ; Szasz2018 are found close to the MIT. It is also distinguished from half-filled $SU(2N\geq 4)$ case, where a direct transition from semimetal to valence bond solid was identifiedZhou2016 ; Li2019 .

Model and Method.—We consider a $SU(4)$ symmetric two-orbital Hubbard model on the honeycomb lattice. The model is given by

[TABLE]

where $\alpha=1,2$ denote two orbits. $c_{i\sigma,\alpha}^{\dagger}$ ($c_{i\sigma,\alpha}$ ) represents the electron creation (annihilation) in orbit $\alpha$ at the $i$th site with spin $\sigma$ ($\sigma=\uparrow,\downarrow$). $n_{i\sigma,\alpha}$ is the electron number operator. For each orbit, the honeycomb lattice on the torus or cylinder is spanned by length vectors $\mathbf{L_{x}}=L_{x}\mathbf{e_{x}}$ and $\mathbf{L_{y}}=L_{y}\mathbf{e_{y}}$, where $\mathbf{e_{x}}=(1,0)$ and $\mathbf{e_{y}}=(1/2,\sqrt{3}/2)$ are two primitive vectors, then the total number of sites for each orbit is $N_{0}=L_{x}\times L_{y}\times 2$ , where $L_{x}$ and $L_{y}$ represent the number of unit cells along the two primitive-vector directions. For the two-orbital system on honeycomb lattice we have the number of sites $N=2\times N_{0}$ (where $2$ denotes two orbits) and the number of electrons $N_{e}=N_{0}$ for a quarter filling.

We set $t$ as the unit of energy and consider Coulomb repulsion $U>0$. Given the additional sublattice and orbital degrees of freedom, we mainly focus on the cylinders or torus with circumferences up to $L_{y}=4$ unit cells ($2\times L_{y}=8$ lattice sites in each orbit). In the present calculations, we keep up to $8000$ states with enough number of sweeps to get the converged data, the truncation error is of the order or less than $10^{-4}$. We also benchmark against quantum Monte Carlo (QMC) for single orbital modelSM .

Metal-Insulator Transition.—We begin with studying the expectation value of the double occupancy $n_{d}$ $n_{d}=\frac{1}{{{N_{0}}}}\sum_{i=1}^{{N_{0}}}{\left\langle{{{\left[{\sum_{\sigma,\alpha}{{n_{i\sigma,\alpha}}}-1}\right]}^{2}}}\right\rangle}$, which corresponds to the first order derivative of the ground-state energy and is a good measurement for the Mott transition. Figures 1 (a) and (b) show $n_{d}$ versus $U/t$ for the model Hamiltonian (1) on torus and cylinder. With the increase of $U/t$, $n_{d}$ monotonically decreases, which eventually freezes the charge degrees of freedom. The kinks of each curve and crossing of different curves with different lattice sizes indicate the transition takes place near $U_{c}/t=2.5\sim 3$, which are independent of lattice geometries we employed. We also show $U\cdot n_{d}$ as a function of $U/t$ in the insets of Fig. 1, where the peaks in the curves may characterize the transition since the derivatives of the curve show discontinuity from both sides of the maximum point.

To establish the nature of each phase, we examine the momentum distribution function $n(\mathbf{k})$, which is defined by the Fourier transformation of single particle propagator $\left\langle{c_{i\sigma}^{\dagger}{c_{j\sigma}}}\right\rangle$, i.e., $n\left(\mathbf{k}\right)=\frac{1}{{{N_{0}}}}\sum\limits_{i,j,\sigma}{\left\langle{c_{i\sigma}^{\dagger}{c_{j\sigma}}}\right\rangle{e^{i{\mathbf{k}}({{\bf{r}}_{\bf{i}}}-{{\bf{r}}_{\bf{j}}})}}}$ . As shown in Figure 2 (a) and (b), $n(\mathbf{k})$ shows distinct behavior before and after the transition. For the metallic phase at $U<U_{c}$, there is a sudden drop near the Fermi surface [see Fig. 2 (a)], while $n(\mathbf{k})$ is near flat without any singularity in the insulating side at $U>U_{c}$ [see Fig. 2 (b)]. To show the evolution of the Fermi surface with increasing $U/t$, we cut the $n(\mathbf{k})$ along the line crossing $\Gamma=(0,0)$ point, as shown in Fig. 2 (c), the smooth change of the line shape indicates the continuity of this MIT, which is consistent with the behavior of $n_{d}$ [see Fig 1]. In addition, We also confirm that results of $n(\mathbf{k})$ are robust against system size as illustrated in Fig. 2 (d).

Density Wave Fluctuations/Orders.— We have identified that the MIT occurs near $U_{c}/t=2.5\sim 3$. Below we examine the electronic and magnetic fluctuations and explore possible orders particularly near the MIT and in the insulating phase. Numerically the most direct evidence to explore the fluctuations or orders is to calculate the correlations and their structure factors. Here, we consider the correlations within the same sublattice, which are the same for the two sublattices. Although the systems on cylinders break translational symmetry due to the open boundaries, we find our results are robust by neglecting the boundary effects when the system size is relatively large.

We first measure the structure factor of charge density-density correlations $N_{\mathbf{q}}=\frac{1}{N_{0}}\sum\limits_{i,j}{\left\langle{n_{i}n_{j}}\right\rangle e^{i\mathbf{q}(\mathbf{r_{i}}-\mathbf{r_{j}})}}$. Figures  3 shows $N_{q}$ in the metallic phase, insulating phase as well as around the vicinity of the MIT. The black dots indicate the data points in the contour plot for finite size system with $L_{x}=12$ and $L_{y}=3$. Here we have substracted the peak at $q=0$ , which is trivially induced by the uniform charge background. Fig. 3 (a) and (d) show that there are no significant peaks in $N_{\mathbf{q}}$ deep inside both the metallic phase at $U<U_{C}$ and insulating phase at $U>U_{C}$, while near the transition point, $N_{\mathbf{q}}$ displays strong charge density wave (CDW) fluctuations, as shown in Fig. 3 (c). When $U/t=3$, the peaks in $N_{q}$ locate at available momenta nearby $\mathbf{M}$ points, implying the incommensurate nature of such CDW fluctuations. We also check the finite size effect by studying wider cylinder with $L_{y}=4$, as shown in the Fig. 5 (a). $N_{q}$ also displays significant peaks around $\mathbf{M}$ but with different wave vectors, which is consistent with the fact that the wave vectors of the incommensurate CDW fluctuations depend on the system geometry. Meanwhile, we also notice that the intensity of $N_{q}$ is weakened with increasing the width of cylinders, as shown in Fig. 3 (c) and Fig. 5 (a), which suggests the CDW fluctuations are not strong enough to form CDW order towards 2D limit.

In the spin channel, we study the static spin structure factor to detect the magnetic fluctuations, which is defined as the Fourier transformation of spin-spin correlations, i.e., $S_{\mathbf{q}}=\frac{1}{N_{0}}\sum\limits_{i,j}{\left\langle{S_{i}^{z}S_{j}^{z}}\right\rangle e^{i\mathbf{q}(\mathbf{r_{i}}-\mathbf{r_{j}})}}$, where $S^{a}_{i}\equiv\sum_{\alpha}\sum_{\sigma\sigma^{\prime}}c^{\dagger}_{i\sigma,\alpha}s^{a}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime},\alpha}$ is the spin operator on site $i$. Due to the $SU(4)$ symmetry of the Hubbard model, the orbital and spin correlations are identical. Figures  4 show the spin/orbital structure factor $S_{\mathbf{q}}$ for different values of $U/t$ for the same $L_{y}=3$ systems with $L_{x}=12$. Both the metallic phase at $U<U_{C}$ [see Fig. 4 (a)] and the insulating phase at $U>U_{C}$ [see Fig. 4 (d)] exhibit the absence of magnetic fluctuations or orders. While near $U_{c}$, $S_{\mathbf{q}}$ displays strong SODW fluctuations, as shown in Fig. 4 (b) and Fig. 4 (c). Interestingly, the SODW fluctuations display competing wave vectors at commensurate momentum $\mathbf{K}$ or $\mathbf{M}$, which depend on the values of ratio between Coulomb interaction and bandwidth. When $U/t=2.5$ [see Fig. 4(c)], there are peaks located at $\mathbf{K^{+}}=(4\pi/3,0)$ and $\mathbf{K^{-}}=(2\pi/3,2\pi/\sqrt{3})$ points in the Brillouin zone. However, further increasing the ratio to $U/t=3$ [see Fig. 4(d)], the peaks of $S(q)$ near $\mathbf{K^{\pm}}$ become broader while much sharper peaks arise at $M=(\pm\pi,\mp\pi/\sqrt{3})$ points, indicating the SODW fluctuations with wave vector $\mathbf{M}$ are more dominant than the ones with wave vector $\mathbf{K}$ near $U_{c}$ This can be seen more clearly by increasing the cylinder width towards 2D, as shown in the Fig. 5 (b) for $S_{q}$ on $L_{y}=4$ cylinders, where the SODW fluctuations display significant peaks around $\mathbf{M}$. In addition, the intensity of $S_{q}$ is also enhanced with increasing the system size [see Fig. 4 (c) and Fig. 5 (b)], indicating the possible existence of SODW order in 2D.

Orbital Polarizations.—We have identified strong charge and spin/orbital dentity wave fluctuations or orders near the quantum phase transition, while the structure factors indicate the featureless nonmagnetic insulating state at $U>U_{c}$ [see Fig. 3(d) and Fig. 4(d) ]. Below, we examine polarization in the spin/orbital channel. To check the possibility of orbital polarizations of the Hamiltonian [see Eq. (1)] at a quarter filling, we define $P=|N_{1}-N_{2}|/N_{0}$, where $N_{\alpha}$ ($\alpha=1,2$) represent the total electron number in orbit $\alpha$, and then we compare the ground-state energies of the unpolarized sector with equal number of electrons in two orbits ($P=0$) and the orbital polarized sector with all electrons in one orbit ($P=1$). As shown in the Fig. 6 for different $L_{y}$ systems, $P=0$ sector always has lower energy than $P=1$ sector, which is independent of the system sizes, implying the absence of the fully orbital polarization for the ground states. Due to the $SU(4)$ symmetry, this result also implies that the ground state is not fully spin-polarized. This conclusion is consistent with the absence of enhancement in the $\mathbf{q}=0$ spin structure factor.

Discussion and Summary.—We now discuss the relevance of our work to real materials. Recently, the two-orbital $SU(4)$ Hubbard model on the honeycomb lattice at quarter filling has been proposed as a realistic model for $\alpha$-ZrCl3 and metal organic frameworks with tricoordinated lattices Yamada . It may be possible to tune the ratio $U/t$ by pressure and to search for spin/orbital density wave and Mott insulator states.

For TBG, the microscopic two-orbital Hamiltonian derived from band structure calculations and Coulomb interactions projected to low-energy bands has more terms than the simplest two-orbital Hubbard model studied here. It includes single-particle hoppings beyond nearest neighborsSenthil , extended density and exchange interactions beyond on-site $U$, pair hoppings Vafek ; Guinea , and the effect of other Moiré bands may be important Balents ; Bernevig ; Vishwanath . These additional terms can stabilize SODWliang-newer or lead to new phases beyond the density wave and Mott insulator states. For example, Ref.Vafek2 found a orbital-polarized state at quarter filling induced by pair hopping, when kinetic energy is set to zero. The fully understanding the physics in TBG requires a more comprehensive study of the microscopic Hamiltonian. Nonetheless, the current two-orbital Hubbard model with the nearest-neighbor hopping and on-site repulsion may be served as a starting and reference point to study correlated electron physics of TBG.

In summary, based on DMRG simulations of the two-orbital Hubbard model on honeycomb lattice, we identify a MIT near $U_{c}/t\approx 2.5\sim 3$ accompanied by commensurate SODW orders and incommensurate CDW fluctuations, and find a nonmagnetic Mott insulator without orbital polarization at $U>U_{c}$. Our results enrich the Mott physics near MIT, where a SODW emerges as an intermediate phase. We also find the SODW is robust against lightly doping and the weak nearest neighbor (NN) interactionsSM ; Wehling2011 . In addition, we have also checked that there is no indication of any density waves at 1/8 fillingSM , which corresponds to one electron per unit cell. Our findings can be tested experimentally in TBG, where applying pressure can drive the system across MIT and deep into the Mott insulator Yankowitz2018 .