Phase diagram of the spin-$\frac{1}{2}$ $J_1$-$J_2$ Heisenberg model on a honeycomb lattice

We use the density matrix renormalization group (DMRG) algorithm to study the phase diagram of the spin-$\frac{1}{2}$ Heisenberg model on a honeycomb lattice with first ($J_1$) and second ($J_2$) neighbor antiferromagnetic interactions, where a $Z_2$ spin liquid region has been proposed. By implementing SU(2) symmetry in the DMRG code, we are able to obtain accurate results for long cylinders with a width slightly over 15 lattice spacings and a torus up to the size $N=2\times 6\times 6$. With increasing $J_2$, we find a Néel phase with a vanishing spin gap and a plaquette valence-bond (PVB) phase with a nonzero spin gap. By extrapolating the square of the staggered magnetic moment $m_s^2$ on finite-size cylinders to the thermodynamic limit, we find the Néel order vanishing at $J_2/J_1\simeq 0.22$. For $0.25<J_2/J_1\le 0.35$, we find a possible PVB order, which shows a fast growing PVB decay length with increasing system width. For $0.22<J_2/J_1<0.25$, both spin and dimer orders vanish in the thermodynamic limit, which is consistent with a possible spin liquid phase. We present calculations of the topological entanglement entropy, compare the DMRG results with the variational Monte Carlo, and discuss possible scenarios in the thermodynamic limit for this region.

Figure 1

Phase diagram of the spin-$\frac{1}{2}$ $J_1$-$J_2$ honeycomb Heisenberg model for $J_2\le 0.35$ obtained by our SU(2) DMRG studies. With increasing $J_2$, the model has a Néel phase for $J_2\lesssim 0.22$ and a PVB phase for $0.25\lesssim J_2\lesssim 0.35$. Between these two phases, there is a small region that exhibits no order in our calculations. The main panel shows Néel order parameter $m_s$ and spin gap $\Delta E_T$. The inset is the sketch of the $J_1$-$J_2$ honeycomb lattice on a $N=2\times L_1\times L_2$ torus (here with four unit cells, $L_1=L_2=4$, along the two primitive vector directions).

Figure 2

Cylinders used in DMRG calculations. (a) ZC4-12 cylinder with zigzag open edges. It has four unit cells along the zigzag direction ($W_y=4\sqrt{3}$) and 12 columns along the axis direction. (b) AC4-12 cylinder with armchair open edges. It has four vertical bonds along the armchair direction ($W_y=6$) and 12 armchair columns along the axis direction. (c) Trimmed ZC cylinder tZC6-18 with trimmed zigzag edges. It has six unit cells along the zigzag direction ($W_y=6\sqrt{3}$) and 18 columns along the axis direction.

Figure 3

Circumference dependence of the ground-state energy per site for $J_2=0.25$ and $0.3$ on torus, AC, tZC, and ZC cylinders. Data are obtained by keeping the optimal states to our computation limit. The truncation errors are below $1\times {10}^{-6}$ except for the largest size at $J_2=0.25$ (tZC9 cylinder with $W_y\simeq 15.588$), where the error is about $5\times {10}^{-6}$. The dashed lines indicate the extrapolations of the energies, which give $-0.4378$ and $-0.4255$ for $J_2=0.25$ and $0.3$, respectively.

Figure 4

(a) $m_s^2$ plotted vs $1/\sqrt{N}$ for the torus clusters $N=2\times 4\times 4$, $2\times 5\times 4$, $2\times 5\times 5$, $2\times 6\times 5$, and $2\times 6\times 6$. The ED data is from Ref. 40. (b), (c) Size dependence of same-sublattice spin structure factor obtained on the torus for $J_2=0.25$ and $0.3$, respectively. The system sizes are $N=2\times 4\times 4$, $2\times 5\times 5$, and $2\times 6\times 6$. (d) $m_s^2$ plotted vs $1/L$ for the ZC$L$-$2L$ cylinder with $L=4,5,6,7,8,9$. Here $m_s^2$ is obtained from $N/2$ spins in the middle part of the sample.

Figure 5

Structure factor of dimer-dimer correlation functions on the $N=2\times 6\times 6$ torus for (a) $J_2=0.25$ and (b) $J_2=0.3$.

Figure 6

The PVB bond texture $B_{i,j}$ for $J_2=0.25$ on an AC6-12 cylinder. The red bonds with negative values have lower NN bond energies. Thus, the red hexagons with negative bonds could indicate the “resonating plaquettes.”

Figure 7

(a) Log-linear plot of ${\epsilon }_{i,j}{B}_{i,j}$ vs distance from the boundary to the bulk for $J_2=0.25$ on AC4-12, AC6-18, and AC8-24 cylinders. The decay length ${\xi }_P$ is obtained by fitting the decay behavior of ${\epsilon }_{i,j}{B}_{i,j}$. (b) Circumference dependence of the decay length ${\xi }_P$ on AC4-12, AC6-18, and AC8-24 cylinders for various $J_2$ couplings.

Figure 8

The PVB bond texture $B_{i,j}$ for $J_2=0.25$ on an AC10-30 cylinder. For clarity, only the left half of the lattice is shown. For this large size, the truncation error is about $5\times {10}^{-6}$ and the bond energies have a some small uncertainty of $\pm 0.001$. At this accuracy, the PVB order vanishes in the middle of the sample.

Figure 9

The PVB bond texture $B_{i,j}$ for $J_2=0.27$ on a tZC6-18 lattice (left half of the lattice is shown). The summation of the bond textures on a hexagon with six negative bond textures (red bonds) is denoted as $E_6$, while that on a hexagon with three red bonds is $E_3$. We define the PVB order parameter as the energy difference between two neighboring such hexagons, $P\equiv |E_6-E_3|$.

Figure 10

(a) Log-linear plot of $P\left(d\right)$ on the tZC9-30 lattice for various $J_2$ couplings. (b) Circumference dependence of decay length ${\xi }_P$ on the tZC3-12, tZC6-18, and tZC9-30 cylinders for various $J_2$ couplings. For $J_2=0.27$ and $0.3$ on the tZC9-30 cylinder, the systems have a long-range PVB order and thus ${\xi }_P$ is divergent. (c) Real-space decay of $P\left(d\right)$ for $J_2=0.3$ on the tZC9-30 cylinder showing bulk PVB order.

Figure 11

The PVB bond texture $B_{i,j}$ for $J_2=0.3$ on the ZC6-24 cylinder with pinning appropriate for the PVB order. We show only the left half part of the lattice. The blue dashed lines indicate the bonds with pinning coupling $J_{\mathrm{pin}}=0.5$.

Figure 12

Log-linear plot of PVB order on the ZC cylinder with PVB pinning $J_{\mathrm{pin}}=0.5$ as shown in Fig. 11 for $J_2=0.25$ and $0.3$.

Figure 13

Spin gap obtained from the torus. (a) At $J_2=0.1$ and $0.15$, the spin gaps are extrapolated to zero as $\Delta E_{T,N}=\alpha /N-\beta /N^{3/2}$ from the samples $N=2\times 3\times 3$, $2\times 4\times 3$, $2\times 4\times 4$, $2\times 5\times 4$, $2\times 5\times 5$, and $2\times 6\times 5$. (b) At $J_2=0.25$ and $0.3$, the spin gaps are extrapolated to finite values as $\Delta E_{T,N}=\Delta E_{T,\infty }+\alpha /N+\beta /N^2$ from the larger samples $N=2\times 4\times 4$, $2\times 5\times 4$, $2\times 5\times 5$, $2\times 6\times 5$, and $2\times 6\times 6$.

Figure 14

Circumference dependence of the entanglement entropy in the large $L_1$ limit measured on both AC and tZC cylinders for (a) $J_2=0.25$ and (b) $J_2=0.3$. The linear extrapolations of the entanglement entropy using data on both cylinders lead to the topological entanglement entropy $\gamma =0.51$ and $0.66$, respectively.

Figure 15

Comparisons of DMRG and VMC results on a torus for (a) spin and (b) dimer correlations at $J_2=0.25$ and $N=2\times 6\times 6$. The VMC wave function is the SPS state of Ref. 44 with $\Delta =0.125$ and $\theta =0.6$ and represents a $Z_2$ spin liquid state (VMC results look similar for a range of $\theta$ including $\theta =0$). Sites $j$ and bonds $(k,l)$ are ordered in a typewriter fashion going first in the ${\vec{a}}_2$ direction in Fig. 1.

Figure 16

Variational energies on the $2\times 5\times 5$ torus comparing the ZF Schwinger boson and SPS slave-fermion Ansätze, together with the exact DMRG ground-state energy. The Dirac spin liquid is obtained by setting $\Delta =0$ in the SPS state. The $U$(1) RVB state is obtained by setting $A_2=0$ in the ZF state. The ZF state can also realize a long-ranged RVB state and can thus provide a good approximation to the Néel state for $J_2<0.2$. In the putative spin liquid region for larger $J_2$, the optimal ZF state is a short-ranged RVB state whose energetics is very similar to the optimal SPS case. Note that the optimal parameters in both the ZF and SPS cases are close to the $U$(1) regime in the respective Ansätze as explained in the text. The DMRG on larger clusters is indispensible in determining the ultimate nature of the ground state.