Phase diagram of interacting spinless fermions on the honeycomb lattice: A comprehensive exact diagonalization study

We investigate the phase diagram of spinless fermions with nearest- and next-nearest-neighbor density-density interactions on the honeycomb lattice at half-filling. Using exact diagonalization techniques of the full Hamiltonian and constrained subspaces, combined with a careful choice of finite-size clusters, we determine the different charge orderings that occur for large interactions. In this regime, we find a two-sublattice Néel-like state, a charge modulated state with a tripling of the unit cell, a zigzag phase, and a charge ordered state with a 12-site unit cell we call Néel domain wall crystal, as well as a region of phase separation for attractive interactions. A sizable region of the phase diagram is classically degenerate, but it remains unclear whether an order-by-disorder mechanism will lift the degeneracy. For intermediate repulsion, we find evidence for a Kekulé or plaquette bond-order wave phase. We also investigate the possibility of a spontaneous Chern insulator phase (dubbed topological Mott insulator), as previously put forward by several mean-field studies. Although we are unable to detect convincing evidence for this phase based on energy spectra and order parameters, we find an enhancement of current-current correlations with the expected spatial structure compared to the noninteracting situation. While for the studied $t-V_1-V_2$ model, the phase transition to the putative topological Mott insulator is preempted by the phase transitions to the various ordered states, our findings might hint at the possibility for a topological Mott insulator in an enlarged Hamiltonian parameter space, where the competing phases are suppressed.

Figure 1

(a) Illustration of the honeycomb lattice with $V_1,V_2$ interactions and the hopping $t$. (b) First (solid line) and second (dashed line) Brillouin zone of the honeycomb lattice including the location of a few special points in the Brillouin zone.

Figure 2

GS energy per site $e_0$ vs $V_2/t$ for (a) $V_1=0$ and (b) $V_1/t=4$ on various clusters. We refer to Appendix pp1 for further details on these finite-size clusters. Lines are guide to the eyes.

Figure 3

Phase diagram in the $(V_1/t,V_2/t)$ parameter space obtained from several exact diagonalization techniques (see text). Dashed lines represent the classical transition lines, see Fig. 4. The semimetal, which is the ground state for noninteracting spinless fermions, has a finite extension in the phase diagram because of its vanishing density of states at the Fermi level. We will argue in the remainder of this article that several other phases can be stabilized for intermediate and/or large interactions: Néel CDW, plaquette/Kekulé (P-K), Néel domain wall crystal (NDWC), zigzag (ZZ) phase, and charge modulation (CM). The region (ST*) is degenerate at the semiclassical level, and it is presently unclear whether and how an order-by-disorder mechanism will lift the degeneracy. Note also the large region of phase separation mostly in the attractive quadrant. Filled symbols correspond to numerical evidence (using level spectroscopy or measurements of correlations, see Sec. 4) obtained mostly on a $N=24$ cluster, which contains the most important points in its Brillouin zone and features the full lattice point group symmetry of the honeycomb lattice. Star symbols denote likely first-order transitions, witnessed by level crossings on the same cluster. Our numerical results do not support any region of topological QAH phase.

Figure 4

Classical phase diagram (a) and a qualitative first order in $t/(V_1,V_2)$ phase diagram (b). The “Zigzag*” and the “Stripy*” phases feature a nontrivial ground-state degeneracy (hence the “*” suffix in “Zigzag*” and “Stripy*”). The three points $V_2=1,V_1=\pm 4$ and $V_2=1,V_1=0$ feature an extensive ground-state degeneracy at the classical level. Upon including the first-order correction due to the finite hopping $t$, two of these points spawn new phases. The $V_2=1,V_1=0$ point develops into the charge modulation (CM) phase, while the $V_2=1,V_1=+4$ point broadens into a novel “Néel domain wall crystal” (NDWC), which is sketched in Fig. 12. All regions and lines in (b) beyond the CM and NDWC phases have no first-order (in $t)$ quantum corrections. Note that the radial direction in the two panels has no physical relevance, i.e., it is not parametrizing $t/V_{1,2}$.

Figure 5

Sketch of a Hamiltonian unit, which renders the classical $V_1-V_2$ problem “frustration free.”

Figure 6

Splitting of the ground-state degeneracy at $V_1/t=-8,V_2/t=6$, yielding the two or six zigzag states as the ground states. Black (red) symbols denote the two different particle-hole symmetry sectors. All levels are singly degenerate unless the (spatial symmetry) degeneracy is indicated by a blue number. The perturbatively generated diagonal terms split the ground space into symmetry related families, whereby the two or six zigzag states become lowest in energy. Note the very small energy difference, indicating a high-order perturbative effect.

Figure 7

Illustration of the stripy charge ordering pattern and the line defects for the $N=24$ (left) and $N=32$ (right) clusters. The left part in each panel indicates a pristine stripy configuration, while the right part shows how another ground state can be reached by a one-dimensional collective transposition of particles. For the $N=24$ cluster, the line move wraps completely around the sample and connects two ground states without defect lines, while for the $N=32$ the indicated shifts of particles generate a ground state with a defect line, where the region along the lines resembles a rotated pristine stripy state. For the $N=24$ cluster, this move is generated in sixth-order perturbation theory, while the move on the $N=32$ cluster is effective in fourth order.

Figure 8

Splitting of the classical ground-state degeneracy at $V_1/t=8,V_2/t=6$ for clusters with $N=12$ (left) and $N=24,24\mathrm{b}$, 28, and 32 sites (right). Black (red) symbols denote the two different particle-hole symmetry sectors. All levels are singly degenerate unless the (spatial symmetry) degeneracy is indicated by a blue number. The order at which processes of the type described in Fig. 7 happen on the clusters are indicated at the bottom of the frame for each system size. In contrast to the case with attractive $V_1$, there is no obvious system-size independent order-by-disorder mechanism at work on the clusters considered.

Figure 9

Low-lying energy spectrum in units of the hopping $t$ for the CM phase at $V_1=0,V_2=1,t/V_2\ll 1$. Data for system sizes $N=18,24,42,54$ are shown. The full arrow denotes the bandwidth of the lowest 18 levels. Note the collapse of the bandwidth for $N=42$ and 54, while the (particle-hole) gap above the 18 states (dashed arrows) remains of order $t$, as illustrated in the inset.

Figure 10

Charge structure factor of the effective model at $V_1=0,V_2=1,t/V_2\ll 1$ plotted in the enlarged Brillouin zone. The filled symbols are data for the ground state at $N=24$ while the empty symbols are for $N=42$. The dashed symbols denote the structure factor in an equal superposition of the $18uud\text{−}ddu$ states for $N=24$. Note the presence of small peaks at momentum ${\mathrm{\Gamma }}^{*}$, indicating a possibly weak sublattice imbalance.

Figure 11

Energetics of the first-order effective model (first order in $t/V_{1,2})$ at $V_1/V_2=4$. (Top) Energy per site in units of $t$ for various cluster sizes. The clusters with 24 and 36 sites have a particularly low energy per site. (Bottom) Single-particle charge gap ${\mathrm{\Delta }}_1\mathrm{p}/t$. This gap amounts to $4\sim 5$ $t$ for the system sizes with a low energy per site.

Figure 12

The Néel domain wall crystal (NDWC): sketch of a classical ground state at $V_1=4,V_2=1$ on the $N=24$ sample, which is maximally flippable within the classical ground-state manifold with respect to the hopping $t$. The shaded regions denote the two Néel domains, and the orange circled bonds along the domain walls are flippable to first order in $t$. The green box indicates a twelve-site unit cell.

Figure 13

Low-lying energy spectrum in units of the hopping $t$ for $V_1=4,V_2=1,t/V_2\ll 1$. Data for system sizes $N=24,36,72$ are shown. The full arrow denotes the bandwidth of the lowest 18 levels $(N=36$ does not have rotational invariance, in that case we only expect six possible ground states). Note the collapse of the bandwidth for larger clusters, while the (particle-hole) gap above the 18 states (dashed arrows) remains approximately constant, as illustrated in the inset.

Figure 14

Charge structure factor of the effective model at $V_1=4,V_2=1,t/V_2\ll 1$ plotted in the enlarged Brillouin zone. The filled (empty) symbols are data for the ground state at $N=24$ $\left(N=72\right)$. The Bragg peaks appear at the $X^{*}$ point, cf. Fig. 1.

Figure 15

Full spectrum vs $V_2/t$ at fixed $V_1/t=0$ labeled with symmetry sectors on $N=24$ cluster. More precisely, states with momentum $\mathrm{\Gamma }$ can be labeled with $C_{6v}$ standard point-group notations; with momentum $X$, only reflection can be used to detect even (e) vs odd (o) states; at the $K$ point, we use $C_{3v}$ point-group notations; at the $M$ point, we can use two reflections to label states that are even/even (e1), even/odd (e2), odd/even (o1), or odd/odd (o2). Open (filled) symbols denote, respectively, states that are even (odd) with respect to particle-hole symmetry.

Figure 16

Same as Fig. 15 at fixed $V_1/t=4$ on the $N=24$ cluster. For intermediate values $V_2/t\simeq 1.5$, low-energy states (indicated with an ellipse) are compatible with a plaquette or Kekulé phase.

Figure 17

Connected density correlation functions between a reference site (open circle) and other sites for $V_1=1$ and $V_1=2$ (at fixed $V_2=0)$ on $N=32$ cluster. Periodic boundary conditions are used. Blue and red correspond to positive/negative correlations.

Figure 18

Same as Fig. 15 as a function of $V_2/t$ at fixed $V_1/V_2=4$ on $N=24$ cluster. The orange boxes denote the two parameter values for which the correlations shown in Fig. 19 were calculated.

Figure 19

Kinetic energy and density correlations computed on the $N=24$ cluster along the $V_1/V_2=4$ line, i.e., for $(V_1/t,V_2/t)$, respectively, equal to $(4,1)$ and $(40,10)$.

Figure 20

Connected density correlations computed on $N=24$ (left) and $N=32$ (right) cluster for $V_1/t=-8$ and $V_2/t=6$. Blue and red correspond to positive/negative correlations.

Figure 21

Kinetic energy correlations at fixed $V_1/t=8$ and $V_2/t=2$ on the $N=42$ cluster. Periodic boundary conditions are used. Blue and red correspond to positive/negative correlations. The width is proportional to the correlation.

Figure 22

Current-current correlations on a given sublattice between a reference bond (black) and other bonds for $V_2/t=1$, 2, and 3 (from left to right) with fixed $V_1=0$ on $N=24$ sample. Periodic boundary conditions are used. Blue and red correspond to positive/negative correlations according to the orientations expected in the QAH phase (taken here to be clockwise on all hexagons). Width is proportional to the correlation.

Figure 23

Same as Fig. 22 for $V_1=0$ and $V_2/t=1$ on various clusters having $C_6$ symmetry. From left to right: $N=26$, 32, 38, and 42.

Figure 24

Same as Fig. 23 for $V_1=0$ and $V_2/t=2$.

Figure 25

Current structure factor ${\mathcal{J}}_q$ as a function of $V_2/t$ for several lattices. Inset: scaling of this structure factor divided by the number of bonds $N_b=3N/2-11$ vs $1/N$ for several $V_2/t$. Open symbols are used for clusters without the $K$ point. Hatched symbols are used for the clusters $N=24$ and 42, which have the $K$ point as well as $C_6$ symmetry at least. All data points are compatible with a vanishing signal in the thermodynamic limit.