Competing states in the SU(3) Heisenberg model on the honeycomb lattice: Plaquette valence-bond crystal versus dimerized color-ordered state

Conflicting predictions have been made for the ground state of the SU(3) Heisenberg model on the honeycomb lattice: Tensor network simulations found a plaquette order [Zhao et al., Phys. Rev. B 85, 134416 (2012)], where singlets are formed on hexagons, while linear flavor-wave theory suggested a dimerized, color-ordered state [Lee and Yang, Phys. Rev. B 85, 100402 (2012)]. In this work we show that the former state is the true ground state by a systematic study with infinite projected-entangled pair states (iPEPS), for which the accuracy can be systematically controlled by the so-called bond dimension $D$. Both competing states can be reproduced with iPEPS by using different unit cell sizes. For small $D$ the dimer state has a lower variational energy than the plaquette state; however, for large $D$ it is the latter which becomes energetically favorable. The plaquette formation is also confirmed by exact diagonalizations and variational Monte Carlo studies, according to which both the dimerized and plaquette states are nonchiral flux states.

Figure 1

Schematic illustration of the two competing states of the SU(3) Heisenberg model on the honeycomb lattice. Thick (thin) bonds correspond to low (high) bond energies. (a) Color-ordered, dimerized state obtained with an 18-site unit cell (shaded in gray). (b) Plaquette state with a six-site unit cell which has no color order.

Figure 2

Local order parameters obtained with iPEPS obtained with two different unit cell sizes. The thickness of the bonds is proportional to the square of the energy of the corresponding bond. Each pie chart on a site visualizes the local color density. (a) The color-ordered, dimerized state obtained with a 18-site unit cell for $D=2$. (b) Plaquette state which has no color order: Each color has the same density on each site ($D=16$).

Figure 3

iPEPS simulation results for the variational energy as a function of inverse bond dimension $1/D$ for different unit cell sizes, compared to the exact diagonalization (ED) results and the lowest energy states from variational Monte Carlo (VMC). For bond dimensions $D\le 8$ the color-ordered, dimerized state has the lowest variational energy, for larger $D$ the plaquette state becomes energetically favorable. The dotted lines are a guide to the eye.

Figure 4

iPEPS simulation results of order parameters as a function of inverse bond dimension $1/D$ of the two competing states. Dotted lines are a guide to the eye. (a) Local ordered moment which is suppressed with increasing bond dimension. For the plaquette state it vanishes for $D\ge 10$, indicating absence of color order. (b) Difference in energy between the strongest bond and the weakest bond in the unit cell, which measures the strength of the plaquette order and dimer order in the two competing states, respectively. In both states the order parameter remains finite also in the infinite-$D$ limit, which shows that the order is stable and that translational symmetry is broken.

Figure 5

Low-energy spectra from ED for $N_s=18$ (top figure) and $N_s=24$ (bottom figure). See the main text for explanation of the different symbols.

Figure 6

Real-space bond energy correlation $\langle P_{ij}{P}_{kl}\rangle -\langle P_{ij}\rangle \langle P_{kl}\rangle$ obtained from ED on clusters of (a) $N_s=18$ and (b) $N_s=24$ sites. The value of the correlation is proportional to the plotted bond strength. Positive (negative) correlations are shown as full blue (dashed red) lines. The reference bond is denoted as the bold black bond.

Figure 7

(a) The 000-flux and (b) the $0\pi \pi$-flux Kekulé states, with six sites and three hexagons in the unit cell. (c) The $\pi \pi \pi$-flux and (d) the $\pi 00$-flux states have 12 sites in the primitive unit cell (shown by the dashed rectangle) or 24 sites in the hexagonal unit cell. These states are characterized by two different absolute values of the hopping amplitudes $t_d$ and $t_h$. The hopping amplitudes on the thin black and gray bonds are $t_d$ and $-t_d$, while on the thick dark and light purple bonds the hopping amplitudes are $t_h$ and $-t_h$, respectively. If the signs of all the hoppings are equal, the product of the hoppings around any of the hexagons is positive: this is the 000-flux state shown in (a). In the other cases, the product of the hoppings around some of the plaquettes is negative, corresponding to a $\pi$ flux on these plaquettes. (e) The chiral $\Phi \Phi \Phi$-flux state with $\Phi =2\pi /3$. The red bonds with arrows denote hoppings $t_h{e}^{i2\pi /3}$ with complex amplitudes. Here one can also introduce the Kekulé modulation for the hoppings by changing the sign on the bonds crossing the boundaries of the unit cell ($t_d$ bonds), resulting in a $\Phi {\Phi }^{\prime }{\Phi }^{\prime }$-flux configuration, with ${\Phi }^{\prime }=5\pi /3$. Also ${\Phi }^{\prime }{\Phi }^{\prime }{\Phi }^{\prime }$- and ${\Phi }^{\prime }\Phi \Phi$-flux configurations can be created by introducing complex hoppings to the $\pi \pi \pi$ and $\pi 00$ cases. (f) Brillouin zone of the honeycomb lattice (black hexagon) with the high-symmetry points $\Gamma =(0,0)$, $\mathrm{K}=(2\pi /3\sqrt{3},2\pi /3)$, and $\mathrm{M}=(0,2\pi /3)$. The Brillouin zone of the 000-, $0\pi \pi$-, $\Phi \Phi \Phi$-, and $\Phi {\Phi }^{\prime }{\Phi }^{\prime }$-flux states is shown by the dark red hexagon, with the high-symmetry points ${\mathrm{K}}_0=(0,4\pi /9)$ and ${\mathrm{M}}_0=(2\pi /3\sqrt{3},2\pi /6)$. The dark green hexagon stands for the Brillouin zone of the $\pi \pi \pi$-, $\pi 00$-, ${\Phi }^{\prime }{\Phi }^{\prime }{\Phi }^{\prime }$-, and ${\Phi }^{\prime }\Phi \Phi$-flux states with ${\mathrm{K}}_{\pi }=(0,2\pi /9)$ and ${\mathrm{M}}_{\pi }=(\pi /3\sqrt{3},\pi /6)$. The index in $\mathrm{M}$ and $\mathrm{K}$ refers to the flux of the central hexagon realized by real hopping amplitudes. The nearest-neighbor distance is chosen to be unity.

Figure 8

(a) Nearest-neighbor bond energy as a function of $t_d/t_h$ in the 72-site cluster. As we change the sign of $t_d/t_h$ from negative to positive, we shift between the $0\pi \pi$-flux and 000-flux states (green squares), or between the $\pi 00$-flux and $\pi \pi \pi$-flux states (purple circles). The energy of the chiral $\Phi \Phi \Phi$-flux state with $\Phi =2\pi /3$ is compared to the $t_d=t_h$ case, while the energy of the chiral $\Phi {\Phi }^{\prime }{\Phi }^{\prime }$-flux state is compared to the $t_d=-t_h$ case. The ${\Phi }^{\prime }{\Phi }^{\prime }{\Phi }^{\prime }$-flux and ${\Phi }^{\prime }\Phi \Phi$-flux states have a higher energy (not shown). The inset of (a) shows the energies around $t_d/t_h=-1$ for the 72- and 288-site clusters. The ground state is degenerate for the $\pi \pi \pi$-flux state when $t_d/t_h>1$, and this is the origin of the scattered energy values. (b) The energies of the $d$ and $h$ bonds versus $t_d$. The hexamerization ($\langle P_h\rangle <\langle P_d\rangle$) is more extended for the $\pi 00$-flux state than for for the $0\pi \pi$-flux state. (c) Schematic drawing of the extension of the hexamerized (plaquette) and dimerized phases that can be read off from the bond energies given in (b). The arrows denote the minima of the energies shown in (a).

Figure 9

(a) The band structures of the free-fermion Hamiltonian for the hexamerized $\pi 00$-flux state ($t_d=-t_h$) along the path ${\mathrm{M}}_{\pi }\Gamma {\mathrm{K}}_{\pi }{\mathrm{M}}_{\pi }$ in the Brillouin zone shown in Fig. 7f. The occupied (dark purple) and unoccupied (light purple) bands are separated by a gap. (b)–(d) The bands of the dimerized $0\pi \pi$-flux state for different values of $t_d$ along the path ${\mathrm{M}}_0\Gamma {\mathrm{K}}_0{\mathrm{M}}_0$. For $t_d>1$ the Fermi surface is a Dirac point at ${\mathrm{K}}_0$—the minimal variational energy corresponds to the case when the Fermi sea touches the $\varepsilon /|t_h|=-1$, $\Gamma$ point for $t_d=-t_h$ [(c)]. The green dashed line shows the Fermi energy.

Figure 10

The expectation value of the $\langle P\left(r\right)\rangle =\langle P_{ij}\rangle$ operator, where $r$ is the distance between the sites $i$ and $j$, for the (a) dimerized $0\pi \pi$- and (b) hexamerized $\pi 00$-flux phases, for the 72- and 288-site clusters with $t_d=-t_h$. The size of the symbols is proportional to the number of bonds having that value of $\langle P\left(r\right)\rangle$. For $r\to \infty$ the value of $\langle P\left(r\right)\rangle$ tends to $1/3$ (denoted by the thin black line), corresponding to the expectation value for the exchange between independent spins, which shows the absence of long-range order. Note the negligible size dependence.

Figure 11

The spin structure factor $S\left(\mathbf{k}\right)$ in the $k$ space for the (a) hexamerized $\pi 00$-flux and (b) dimerized $0\pi \pi$-flux states, calculated from the 648-site cluster, with $t_d=-t_h$. The light green hexagons indicate the Brillouin zone of the honeycomb lattice, with the high-symmetry points $\Gamma$, $\mathrm{K}$, and $\mathrm{M}$. The $S\left(\mathbf{k}\right)$ for the $0\pi \pi$-flux state is peaked one-third of the way on the line between the points ${\mathrm{M}}_{▴}$ and ${\mathrm{K}}_{▴}$ in the second Brillouin zone of the honeycomb lattice, as can also be seen in (c), where the $S\left(\mathbf{k}\right)$ is shown along the path $\Gamma {\mathrm{M}}_{▴}{\mathrm{K}}_{▴}\Gamma$ for the 288- and 648-site clusters. The $S\left(\mathbf{k}\right)$ for the $\pi 00$-flux phase (open symbols), also with $t_d=-t_h$, shows no singular behavior. The spin structure factors for the two flux states are nearly identical away from the ${\mathrm{M}}_{▴}$ point.

Figure 12

The stability of the dimerized and hexamerized states versus the formation of long-range order (LRO), in the (a) $0\pi \pi$- and (b) $\pi 00$-flux states, respectively. For both cases the energy is minimal for ${\epsilon }_{\text{LRO}}=0$ and $t_{\text{LRO}}=t_h$, where there is no long-range order. The calculations were made with step sizes $\Delta t_{\text{LRO}}=0.1\left|t_h\right|$ and $\Delta {\epsilon }_{\text{LRO}}=0.2\left|t_h\right|$. For $t_{\text{LRO}}=0$ the fermionic band structure collapsed to a few highly degenerate bands, which caused difficulties in the Monte Carlo code we used for the calculation.