Transitions to valence-bond solid order in a honeycomb lattice antiferromagnet

We use quantum Monte Carlo methods to study the ground-state phase diagram of a $S=1/2$ honeycomb lattice magnet in which a nearest-neighbor antiferromagnetic exchange $J$ (favoring Néel order) competes with two different multispin interaction terms: a six-spin interaction $Q_3$ that favors columnar valence-bond solid (VBS) order, and a four-spin interaction $Q_2$ that favors staggered VBS order. For $Q_3\sim Q_2\gg J$, we establish that the competition between the two different VBS orders stabilizes Néel order in a large swath of the phase diagram even when $J$ is the smallest energy scale in the Hamiltonian. When $Q_3\gg (Q_2,J)$ $\mathrm{[}{Q}_2\gg (Q_3,J)\text{]}$, this model exhibits at zero temperature phase transition from the Néel state to a columnar (staggered) VBS state. We establish that the Néel-columnar VBS transition is continuous for all values of $Q_2$, and that critical properties along the entire phase boundary are well characterized by critical exponents and amplitudes of the noncompact ${\mathrm{CP}}^1$ $\left({\mathrm{NCCP}}^1\right)$ theory of deconfined criticality, similar to what is observed on a square lattice. However, a surprising threefold anisotropy of the phase of the VBS order parameter at criticality, whose presence was recently noted at the $Q_2=0$ deconfined critical point, is seen to persist all along this phase boundary. We use a classical analogy to explore this by studying the critical point of a three-dimensional $XY$ model with a fourfold anisotropy field which is known to be weakly irrelevant at the three-dimensional $XY$ critical point. In this case, we again find that the critical anisotropy appears to saturate to a nonzero value over the range of sizes accessible to our simulations.

Figure 1

(a) Columnar VBS order on the honeycomb lattice: dark links represent higher values of $\langle P_l\rangle$ (the singlet projection operator on this link) than light links. If dark links are instead reinterpreted as representing lower values of $P_l$, one obtains a representation of plaquette VBS order at the same wave vector. (b) Staggered VBS order on the honeycomb lattice, where again dark links represent higher values of $\langle P_l\rangle$. In both panels, we have created different ordered domains (represented by different colors and separated by dashed lines) by introducing a defect. As already discussed [13], the defect has a spinful core (a free spin $1/2$ sits at the domain wall intersection) for columnar/plaquette VBS whereas the core is spinless for the staggered VBS. Also shown are our conventions for labeling unit cells $\vec{r}$, bonds $\mu$ belonging to unit cells $\vec{r}$, and $A$ and $B$ sublattice sites in unit cell $\vec{r}$. (c) Schematic representations of the four- and six-spin interaction terms $Q_2$ and $Q_3$.

Figure 2

Color maps of Binder cumulants (top: Néel Binder cumulant $g_M$, middle: VBS columnar Binder cumulant $g_{\psi }$, bottom: staggered Binder cumulant $g_{\varphi })$ for different values in the $\left(Q_2,Q_3\right)$ parameter space, for a system of linear size $L=24$. Low values indicate long-range order, while high values indicate absence of order. These results allow us to map the phase diagram where Néel, columnar, and staggered VBS phases can be identified (green lines are indications of approximate phase boundaries).

Figure 3

The first-order transition from Néel to staggered VBS order is readily identified by the sharp jumps in the corresponding order parameters (filled squares for $\langle M^2\rangle$ shown on left $y$ axis, filled circles for $\langle {\varphi }^2\rangle$ shown on right $y$ axis), for three values of $Q_3$ (system size $L=24)$.

Figure 4

Crossing plot of Binder cumulants for different system sizes: Néel cumulant $g_M$ (top panel, for $Q_2=0.14)$ and columnar VBS cumulant $g_{\psi }$ (bottom panel, for $Q_2=0.60)$. Symbols are QMC data, solid lines fits to the finite-size scaling form (see text). For the fits, a particular choice of critical window, minimum system size included, and order of universal function has been shown here which gave ${\chi }^2$ per degree of freedom equal to 1.53 and 0.97 for the plots respectively. For estimates on overall error bars, refer to Table 1.

Figure 5

Scaling collapse for different system sizes of the Néel order parameter $\langle M^2\rangle$ (top panel, $Q_2=0.85)$ and columnar VBS order parameter $\langle |\psi {|}^2\rangle$ (bottom panel, $Q_2=20.0)$ in the critical region. Critical point $Q_{3c}$ and critical exponents are obtained by fits to the standard finite-size scaling forms (see text). For the fits, again a particular choice of critical window, minimum system size included, and order of universal function has been shown here which gave ${\chi }^2$ per degree of freedom equal to 1.25 and 1.15 for the plots respectively. For estimates on overall error bars, refer to Table 1.

Figure 6

Finite-size dependence of $W_3$ close to the critical point for $Q_2=0.14$ (top panel), $Q_2=0.60$ (middle panel), $Q_3=20.0$ (bottom panel). In each case, we display data for one value of $Q_3$ closer to our estimate of $Q_{3c}$, one slightly above and one slightly below.

Figure 7

3d XY model with fourfold anisotropic field: crossing plot for different system sizes for the Binder cumulant $B$ (top panel, for $h=0.05)$ and correlation ratio $R$ (bottom panel crossing plot, for $h=2)$.

Figure 8

3d XY model with fourfold anisotropic field: collapse of the order parameter $\langle \left|m\right|\rangle$ (top panel, for $h=0)$ and the Binder cumulant $B$ (bottom panel, for $h=1)$, according to the scaling forms mentioned in the text. For estimates on overall error bars, refer to Table 2.

Figure 9

System-size dependence of the anisotropy parameter $W_4$ close to criticality for three different temperatures: above, below, and very close to the critical temperature $T_c$. Top panel: $h=0.5$, bottom panel: $h=2$.

Figure 10

Scaling collapse plots according to the scaling form $W_3=g_{W_3}\left((Q_3-Q_{3c})L^{1/{\nu }_3}\right)$ for the dimensionless anisotropy quantifier $W_3$ in the $J\text{−}{Q}_2\text{−}{Q}_3$ model for $Q_2=0.14$ (top panel) and $Q_2=20$ (bottom panel). For the fits, similar to Sec. 4a, a particular choice of critical window, minimum system size included, and order of universal function was taken here which gave ${\chi }^2$ per degree of freedom equal to 1.26 and 1.31 for the plots respectively.

Figure 11

Scaling collapse plots according to the scaling form $W_4=g_{W_4}\left((T-T_c)L^{1/{\nu }_4}\right)$ for the dimensionless anisotropy quantifier $W_4$ in the 3d XY model with fourfold anisotropy field for $h=0.5$ (top panel) and $h=2$ (bottom panel). For estimates on overall error bars, refer to Table 4.

Figure 12

3d XY model with threefold anisotropic field: temperature dependence of the Binder cumulant (top panels) and threefold anisotropy quantifier $W_3$ (bottom panels) for two different values of $h_3=0.5,1.0$.

Figure 13

3d XY model with fivefold anisotropic field: temperature dependence of the Binder cumulant (top panels) and fivefold anisotropy quantifier $W_5$ (bottom panels) for two different values of $h_5=1,10$.