Projected wave function study of ${\mathbb{Z}}_2$ spin liquids on the kagome lattice for the spin-$\frac{1}{2}$ quantum Heisenberg antiferromagnet

Motivated by recent density-matrix renormalization group (DMRG) calculations [Yan, Huse, and White, Science 332, 1173 (2011)], which claimed that the ground state of the nearest-neighbor spin-$1/2$ Heisenberg antiferromagnet on the kagome lattice geometry is a fully gapped spin liquid with numerical signatures of ${\mathbb{Z}}_2$ gauge structure, and a further theoretical work [Lu, Ran, and Lee, Phys. Rev. B 83, 224413 (2011)], which gave a classification of all Schwinger-fermion mean-field fully symmetric ${\mathbb{Z}}_2$ spin liquids on the kagome lattice, we have thoroughly studied Gutzwiller-projected fermionic wave functions by using quantum variational Monte Carlo techniques, hence implementing exactly the constraint of one fermion per site. In particular, we investigated the energetics of all ${\mathbb{Z}}_2$ candidates (gapped and gapless) that lie in the neighborhood of the energetically competitive U(1) gapless spin liquids. By using a state-of-the-art optimization method, we were able to conclusively show that the U(1) Dirac state is remarkably stable with respect to all ${\mathbb{Z}}_2$ spin liquids in its neighborhood, and in particular for opening a gap toward the so-called ${\mathbb{Z}}_2[0,\pi ]\beta$ state, which was conjectured to describe the ground state obtained by the DMRG method. Finally, we also considered the addition of a small second nearest-neighbor exchange coupling of both antiferromagnetic and ferromagnetic type, and obtained similar results, namely, a U(1) Dirac spin-liquid ground state.

Figure 1

The ${\mathbb{Z}}_2[0,\pi ]\beta$ SL Ansatz; black (grey) bonds denote 1st NN real hopping (2nd NN real hopping and real spinon pairing) terms; solid (dashed) black bonds have ${\mathrm{s}}_{ij}=1$ ($-1$), solid (dashed) grey bonds have ${\nu }_{ij}=1$ ($-1$), see Eq. (4). The 1st NN (2nd NN) mean-field Ansatz is written as $U_{\langle ij\rangle }=\pm {\sigma }_3$ [$U_{\langle \langle ij\rangle \rangle }=\pm ({\chi }_2{\sigma }_3+{\Delta }_2{\sigma }_1)$]. The SU(2) flux $P$, through elementary triangles (e.g., 123) is $P_{123}={\sigma }_3$, and that through triangles formed by two 1st NN and one 2nd NN bonds (e.g., 234) is $P_{234}=-$(${\chi }_2{\sigma }_3+{\Delta }_2{\sigma }_1$). Their commutator is nonzero, $[P_{123},P_{234}]=(-2i{\sigma }_2){\Delta }_2$. Hence, a finite ${\Delta }_2$ breaks the U(1) gauge structure down to ${\mathbb{Z}}_2$, and opens up an energy gap via the Anderson-Higgs mechanism.[20, 23]

Figure 2

A typical variational Monte Carlo optimization run for the ${\mathbb{Z}}_2[0,\pi ]\beta$ wave function: (a) variational parameters ${\Delta }_2$, ${\chi }_2$, $\mu$, and ${\zeta }_{\mathrm{R}}$ and (b) energy, as a function of SR iterations. In (a), the initialized parameter values are ${\Delta }_2={\chi }_2=1$, $\mu =-0.8$, and ${\zeta }_{\mathrm{R}}=0.3$. The U(1) 2nd NN $[0,\pi ;0,\pi ]$ Dirac SL corresponds to ${\Delta }_2=0$, ${\chi }_2=-0.0186\left(2\right)$, ${\zeta }_{\mathrm{R}}=0$, as found in Ref. 26. The optimized parameter values are obtained by averaging over a much larger number of converged SR iterations than shown above. In (c), the variation in energy upon addition of a small ${\Delta }_2$ (both for ${\zeta }_{\mathrm{R}}=0$, and optimized ${\zeta }_{\mathrm{R}}$ for each value of ${\Delta }_2$) upon the $[0,\pi ;0,\pi ]$ Dirac SL is shown, the increase in energy is apparent.