Resonating valence-bond physics on the honeycomb lattice

We study bond and spin correlations of the nearest-neighbor resonating valence bond (RVB) wave function for a SU(2) symmetric $S=1/2$ antiferromagnet on the honeycomb lattice. We find that spin correlations in this wave function are short ranged, while the bond energy correlation function takes on an oscillatory power-law form $D\left(\vec{r}\right)\sim cos(\mathbf{Q}·\vec{r})/{\left|\vec{r}\right|}^{{\eta }_w\left(2\right)}$, where $\mathbf{Q}=(2\pi /3,-2\pi /3)$ is the wave vector corresponding to “columnar” valence-bond solid order on the honeycomb lattice, and ${\eta }_w\left(2\right)\approx 1.49\left(3\right)$. We use a recently introduced large-$g$ expansion approach to relate bond-energy correlators of the $\mathrm{SU}\left(g\right)$ wave function to dimer correlations of an interacting fully packed dimer model with a three-dimer interaction of strength $V\left(g\right)=-ln(1+1/g^2)$. Putting $g=2$, we find numerically that the dimer correlation function $D^d\left(\vec{r}\right)$ of this dimer model has power-law behavior $D^d\left(\vec{r}\right)\sim cos(\mathbf{Q}·\vec{r})/{\left|\vec{r}\right|}^{{\eta }_d\left(2\right)}$ with ${\eta }_d\left(2\right)\approx 1.520\left(15\right)$, in rather good agreement with the wave function results. We also study the same quantities for $g=3,4,10$ and find that the bond-energy correlations in the $\mathrm{SU}\left(g\right)$ wave function are consistently well reproduced by the corresponding dimer correlations in the interacting dimer model.

Figure 1

$\left(1\right)$ and $\left(2\right)$ represent the two flippable configurations of dimers around a hexagon (here dimers are represented by ellipses). The weight of configuration $\left(1\right)$ in the interacting dimer model gets contributions from the two loop configurations $\left(a\right)$ and $\left(b\right)$ of the loop model equivalent to the nnRVB wave function. This is captured by the effective interactions between dimers worked out in the text.

Figure 2

$D_2$ shows a configuration with two flippable hexagons that share a dimer. The weight for such a configuration in the interacting dimer model gets contributions from loop configurations $L_3$, $L_4$, and $L_5$. Similarly, $D_3$ represents one of the two flippable arrangements of dimers around on the perimeter of a double hexagon. The weight for such a configuration in the interacting dimer model gets contributions from loop configurations $L_6$, $L_7$, and $L_4$.

Figure 3

The honeycomb lattice is constructed using a two point basis of sites $A$ and $B$ “belonging” to each Bravais lattice point $\vec{r}=m{\widehat{e}}_x+n{\widehat{e}}_y$ with integer $m$ and $n$, which, in our convention, is the coordinate of the $B$-sublattice site. When talking of bond energies or dimer occupation numbers, we use the convention that three types of bonds 0, 1, and 2 “belong” to each $\vec{r}$, as shown in the figure.

Figure 4

$B$-sublattice honeycomb sites are labeled by $\vec{r}$ and their corresponding triangular lattice counterparts are labeled $\vec{R}$. Also shown here are the lines across which dimer or valence-bond flux is calculated to get three winding numbers $w_x$, $w_y$, and $w_z$ as described in the text.

Figure 5

Spin correlation function between two sites on the same sublattice, separated by $\vec{r}$. The fit is to an exponentially decaying function $cexp(-r/\xi )$, with the best-fit value of correlation length $\xi =0.550\left(14\right)$.

Figure 6

$D\left({\vec{R}}_{L/6}\right)$ as a function of $L$ for $g=2,3,4$ and 10. ${\eta }_w\left(g\right)$ extracted from a fit to the form $c{L}^{-{\eta }_w\left(g\right)}$ gives ${\eta }_w\left(2\right)=1.49\left(3\right)$, ${\eta }_w\left(3\right)=1.74\left(6\right)$, ${\eta }_w\left(4\right)=1.88\left(5\right)$, and ${\eta }_w\left(10\right)=2.04\left(7\right)$. ${\eta }_w\left(g\right)$ is expected to equal $1/\left({\kappa }_w\left(g\right)\sqrt{3}\right)$, providing us a way of estimating ${\kappa }_w\left(g\right)$.

Figure 7

$N\left(\vec{q}\right)$ is the correlator of the average valence-bond occupation variables (defined in Sec. 3) in the limit of small $|\vec{q}|$, measured in the zero-winding sector (defined in Sec. 3) of our wave function simulations. The extrapolation to $|\vec{q}|\to 0$ yields intercepts of $0.2046\left(2\right)$, $0.2396\left(2\right)$, $0.2546\left(3\right)$, and $0.2722\left(2\right)$ for $g=2,3,4$, and 10, respectively. These intercepts are expected to equal $1/\left(4\pi {\kappa }_w\right)$, and provide us an accurate estimate of ${\kappa }_w\left(g\right)$.

Figure 8

Winding sector probabilities for the nnRVB simulation at $g=2,3,4$ and 10, plotted as a function of $w^2\equiv w_x^2+w_y^2+w_z^2$. The fit is to an exponentially decaying function $aexp[-c\left(g\right)w^2]$, with best-fit values $c\left(2\right)=1.218\left(5\right)$, $c\left(3\right)=1.045\left(1\right)$, $c\left(4\right)=0.982\left(3\right)$, and $c\left(10\right)=0.919\left(3\right)$. $c\left(g\right)$ is expected to equal $\pi {\kappa }_w\left(g\right)$, providing a means of estimating ${\kappa }_w\left(g\right)$.

Figure 9

$D^{\left(d\right)}\left({\vec{R}}_{L/6}\right)$ as a function of $L$ for $g=2,3,4$ and 10. ${\eta }_d\left(g\right)$ extracted from a fit to the form $c{L}^{-{\eta }_d\left(g\right)}$ gives ${\eta }_d\left(2\right)=1.520\left(15\right)$, ${\eta }_d\left(3\right)=1.79\left(2\right)$, ${\eta }_d\left(4\right)=1.85\left(4\right)$, and ${\eta }_d\left(10\right)=1.96\left(9\right)$. ${\eta }_d\left(g\right)$ is expected to equal $1/\left({\kappa }_d\left(g\right)\sqrt{3}\right)$, providing us a way of estimating ${\kappa }_d\left(g\right)$.

Figure 10

$N\left(\vec{q}\right)$ is the correlator of the average dimer correlation function (defined in Sec. 3) in the limit of small $|\vec{q}|$, measured in the zero-winding sector (defined in Sec. 3) of our dimer-model simulations. The extrapolation to $|\vec{q}|\to 0$ yields intercepts of $0.2074\left(1\right)$, $0.2408\left(1\right)$, $0.2551\left(1\right)$, and $0.2721\left(1\right)$ for $g=2,3,4$, and 10, respectively. These intercepts are expected to equal $1/\left(4\pi {\kappa }_d\right)$ and provide us an accurate estimate of ${\kappa }_d\left(g\right)$.

Figure 11

Winding sector probabilities for the dimer-model simulation at $g=2,3,4$ and 10, plotted as a function of $w^2\equiv w_x^2+w_y^2+w_z^2$. The fit is to an exponentially decaying function $aexp[-c\left(g\right)w^2]$, with the best-fit values $c\left(2\right)=1.211\left(4\right)$, $c\left(3\right)=1.036\left(4\right)$, $c\left(4\right)=0.979\left(4\right)$, and $c\left(10\right)=0.918\left(3\right)$. $c\left(g\right)$ is expected to equal $\pi {\kappa }_d\left(g\right)$, providing a means of estimating ${\kappa }_d\left(g\right)$.

Figure 12

Comparisons of values of ${\kappa }_d$ obtained using different methods for $g=2,3,4$, and 10.

Figure 13

Comparisons of values of ${\kappa }_w$ obtained using different methods for $g=2,3,4$, and 10.

Figure 14

Comparisons of the estimated values of ${\kappa }_w$ and ${\kappa }_d$, obtained from our wave function and dimer model studies for $g=2,3,4,10$. Our consolidated estimates, considering all three ways of extracting $\kappa$ on an equal footing, are as follows: ${\kappa }_d\left(2\right)=0.383\left(5\right)$, ${\kappa }_w\left(2\right)=0.388\left(6\right)$, ${\kappa }_d\left(3\right)=0.328\left(8\right)$, ${\kappa }_w\left(3\right)=0.331\left(8\right)$, ${\kappa }_d\left(4\right)=0.312\left(7\right)$, ${\kappa }_w\left(4\right)=0.311\left(8\right)$, ${\kappa }_d\left(10\right)=0.292\left(14\right)$, and ${\kappa }_w\left(10\right)=0.290\left(10\right)$. Error bars in these consolidated estimates reflect the spread between three different ways of estimating ${\kappa }_w$ and ${\kappa }_d$, as well as statistical errors in the individual estimates.