Valence-bond crystals in the kagomé spin-1/2 Heisenberg antiferromagnet: a symmetry classification and projected wave function study

In this paper, we do a complete classification of valence-bond crystals (VBCs) on the kagomé lattice based on general arguments of symmetry only and thus identify many new VBCs for different unit cell sizes. For the spin-1/2 Heisenberg antiferromagnet, we study the relative energetics of competing gapless spin liquids (SLs) and VBC phases within the class of Gutzwiller-projected fermionic wave functions using variational Monte Carlo techniques, hence implementing exactly the constraint of one fermion per site. By using a state-of-the-art optimization method, we conclusively show that the U(1) Dirac SL is remarkably stable towards dimerizing into all 6-, 12- and 36-site unit cell VBCs. This stability is also preserved on addition of a next-nearest-neighbor super-exchange coupling of both antiferromagnetic and ferromagnetic (FM) type. However, we find that a 36-site unit cell VBC is stabilized on addition of a very small next-nearest-neighbor FM super-exchange coupling, i.e. |J2| ≈ 0.045, and this VBC is the same in terms of space-group symmetry as that obtained in an effective quantum dimer model study. It breaks reflection symmetry, has a nontrivial flux pattern and is a strong dimerization of the uniform RVB SL.


Introduction
For many decades, physicists have been actively searching for playgrounds that are 'hot' enough to melt magnetic freezing at temperatures well below the characteristic interaction energy scales in the system. This melting being fueled by quantum fluctuations leads to stabilization of exotic quantum paramagnetic phases of matter [1]. Representatives of such phases are spin liquids (SLs) and valence-bond crystals (VBCs); the former preserve lattice symmetries and the latter break them, according to a generally accepted definition. Long before any experimental hints, theoreticians such as Pomeranchuk already conjectured the existence of SLs [2], which were later advocated by Anderson to be possible appropriate ground states for the spin-1/2 Heisenberg antiferromagnet [3,4]. On the experimental side, the drought in the search for SLs ended with the discovery of Herbertsmithite (ZnCu 3 (OH) 6 Cl 2 ), a compound with perfect kagomé lattice geometry, belonging to the paratacamite family [5,6,7,8,9,10,11,12,13]. In it, the combination of low spin value (S = 1/2), low-dimensionality (d = 2) and coordination number (z = 4) and frustrating nearest-neighbor (NN) antiferromagnetic (AF) super-exchange interactions on a non-bipartite lattice leads to amplification of quantum fluctuations that stabilize a quantum paramagnet. Indeed, all experimental probes on Herbertsmithite point to a SL behavior down to 20 mK (∼ J/10 4 ), which was established on the magnesium version of Herbertsmithite (i.e. MgCu 3 (OH) 6 Cl 2 ) [14,15,16]. Furthermore, Raman spectroscopic studies on Herbertsmithite hint at a gapless (algebraic) SL [17].
In this work, we will study these non-magnetic phases within a Schwinger fermion formulation of the spin model. Within this approach, the projected gapless (algebraic) U(1) Dirac SL has the best variational energy [32]; despite being a marginally stable phase, it was argued in [33] to be stable against a certain class of perturbations. Explicit numerical calculations using projected wave functions have in fact shown it to be stable (locally and globally) w.r.t. dimerizing into all known VBC perturbations [32,34,35]. Furthermore, it was shown that within this class of Gutzwiller projected wave functions, all the fully symmetric gapped Z 2 SLs have a higher energy compared to the U(1) Dirac SL [36]. Similar conclusions were also reached within the Schwinger boson approach to the spin model [56,57]. Note that a simple tensor network (PEPS) representation of such a projected bosonic RVB ansatz can be constructed and has been studied in [58].
In this paper, in section 2 we first perform a systematic symmetry classification of VBC patterns on the kagomé lattice and thus identify and enumerate many new VBCs, independent of the formalism used to study them. In section 3, we address the question of relative energetics of SL and VBC phases. In particular, in section 3.1.1 we show that the U(1) Dirac SL is remarkably stable w.r.t. dimerizing into any of these new VBCs. This stability is also preserved upon addition of a finite next-nearest-neighbor (NNN) super-exchange coupling of both AF and ferromagnetic (FM) type. Such a NNN coupling might be a possible perturbation in Herbertsmithite. In section 3.1.2, we show that a broken symmetry phase is stabilized on addition of a small NNN FM coupling, which is consistent with the findings in [59]. This VBC has a 36-site unit cell with a nontrivial flux pattern threading its plaquettes and it is found to be a strong dimerization of another competing U(1) gapless SL, the so-called uniform RVB SL [32]. This 36-site unit cell VBC has a lower symmetry as compared to that studied in our previous work [35] and has precisely the same symmetry as that identified in QDM studies [49,50,55]. Thus, here we mainly establish the stability of the U(1) Dirac SL w.r.t. an extremely large class of potential VBC instabilities and detect a non-trivial 36-site unit cell VBC instability of the uniform RVB SL which is stabilized on addition of a very weak NNN FM super-exchange coupling to the Hamiltonian.

The model, wave functions and the numerical technique
The Hamiltonian for the spin-1/2 quantum Heisenberg J 1 −J 2 model iŝ where ij and ij denote sums over NN and NNN pairs of sites, respectively. Thê S i are spin-1/2 operators at each site i. In the following, we will consider J 1 > 0 (AF) and both FM and AF super-exchange J 2 ; all energies will be given in units of J 1 . The physical variational wave functions are defined by projecting noncorrelated fermionic states: where P G = i (1 − n i,↑ n i,↓ ) is the full Gutzwiller projector enforcing the one fermion per site constraint. Here, |Ψ MF (χ ij , ∆ ij , µ, ζ) is the ground state of a mean-field  Hamiltonian constructed out of Schwinger fermions and containing hopping, chemical potential and singlet pairing terms: where χ ij = χ * ji and ∆ ij = ∆ ji . Besides the chemical potential µ, we will also consider real and imaginary components of on-site pairing, which are absorbed in ζ.
The SL phases are characterized by different patterns of distribution of underlying SU(2) gauge fluxes through the plaquettes which are implemented by a certain distribution of the phases of χ ij and ∆ ij on the lattice links. Since in a SL state |χ ij | 2 + |∆ ij | 2 is constant for each geometrical distance, a complete specification of a SL state up to nth NN amounts to specifying, in addition to the SU(2) fluxes, the optimized magnitude of hopping and pairing parameters at each geometrical distance and the specification of the on-site terms µ and ζ [60,61]. On the other hand, in a VBC state |χ ij | 2 + |∆ ij | 2 may be different from bond to bond and, therefore, the specification of VBCs amounts to giving the pattern of amplitudes of χ ij and ∆ ij at each geometrical distance, in addition to specifying the SU(2) fluxes through the plaquettes.These parameters are the ansätze of a given state and serve as the variational parameters in the physical wave function that are optimized within the variational Monte Carlo scheme to find the energetically best state. It is worth mentioning that we use a sophisticated implementation of the stochastic reconfiguration (SR) optimization method which allows us to obtain an extremely accurate determination of variational parameters [62,63]. Indeed, small energy differences are effectively computed by using a correlated sampling, which makes it possible to strongly reduce statistical fluctuations. This feature is especially important for the spin-1/2 QHAF since the energies of all the competing phases are rather close.  The square boxes contain the VBC names followed by their respective symmetry PG. The 'parent' (maximally symmetric) VBCs are marked in red and those which have been found as competing ground states in studies using quantum dimer models are marked in pink [49,50,55]. The corresponding VBC patterns, and their discussion, are given in the text. As much as possible, we use labeling for the VBCs which is similar to that used in [55].

Parent spin liquid states
The ansatz for the energetically best variational state, the U(1) Dirac SL, is given in figure 1a. Due to the U(1) flux ϕ being 0 and π [exp (iϕ) = plaquette χ ij ] through triangles and hexagons, respectively, it is denoted as [0, π] SL. In its mean-field band structure the Fermi surface collapses to two points at which the spectrum becomes relativistic with Dirac conical excitations [32]. Another energetically competing state, the uniform RVB SL, has zero flux through all plaquettes and is therefore denoted as [0, 0] SL. Its mean-field band structure consists of large circular spinon Fermi surfaces [34]. Both these states are fully symmetric, U(1) gapless SLs and can be extended to include second NN hoppings, leading to a gain in energy without changing Figure 3: The most symmetric 12-site unit cell VBCs: the center of symmetry is marked 'C' (the center of the shaded hexagon), around which bonds connected by the given PG symmetry operations are marked with the same color and style of the line. We will henceforth refer to these bonds as being in the same class. (a) The Star-VBC has the maximal PG symmetry, C 6v ; hence it acts as a 'parent' VBC. Its bonds breakup into three distinct classes. (b) The VBC 0 lacks crystallographic axes reflection symmetries in contrast to the SVBC; thus its symmetry group is reduced to C 6 . It has four classes of bonds. (c) The Star-VBC α has reduced (2π/3) rotation symmetry but preserves reflection symmetry; thus its symmetry group is C 3v . It has six classes of bonds.
their nature [35]. It is worth noting that the effect of projection on these mean-field states can be drastic.

Symmetry classification and enumeration of valence-bond crystals (VBCs)
The VBC states on the kagomé lattice break its elementary (three-site) unit cell translation symmetry with different unit cell sizes, which describe their modulation.
In previous studies [51,52,53,54,43,44,45,47,49,50,55], using different methods, VBCs with 6-, 12-, 18-and 36-site unit cells were identified as possible ground states of the spin-1/2 QHAF. In this work, we will restrict our analysis to VBCs with 6-, 12-and 36-site unit cells. For each unit cell size with a given center of symmetry, we enumerate VBCs starting from the maximally symmetric (C 6v ) 'parent' VBC and systematically break point group symmetry elements, right down to the VBC with no symmetry at all. This results in an enumeration of 19 VBCs in total, 9 VBCs each for 12-and 36-site unit cells and 1 VBC for the 6-site unit cell (see figure 2). Only 6 out of the 19 VBC have been studied previously. In this paper, we will study the possibility of any of these VBCs to occur as the ground state. Figure 4: Other 12-site unit cell VBCs: (a) the VBC α 0 has only reduced rotation symmetry (2π/3); thus in contrast to SVBC α its symmetry group is reduced to C 3 . It has eight classes of bonds. (b) The diamond-VBC has two perpendicular axes of reflection symmetry, thus giving rise to C 2v symmetry, with seven classes of bonds. (c) The VBC 3 has only π-rotation symmetry; thus its symmetry group is C 2 . It has 12 classes of bonds. (d) The VBC 1 possesses only a single axis of reflection symmetry which bisects the sides of the shaded hexagon; consequently, its symmetry group is C 1v . It has 14 classes of bonds. (e) The VBC 1 has the same symmetry as VBC 1 , but its reflection symmetry axis passes through a vertex of the shaded hexagon; we shall denote the symmetry group as C 1v to distinguish it from that of VBC 1 . It has 12 classes of bonds. (f) The VBC 2 has no symmetry whatsoever; hence its symmetry group is just the identity, denoted here as E. Consequently, it has 24 distinct classes of bonds.

12-site unit cell VBCs
The kagomé lattice can be viewed as a triangular lattice of 12-site blocks shaped in the form of 'stars'. Within this picture, it was argued in [53,54] that the ground state of the spin-1/2 QHAF has possible long-range singlet order that settles in this triangular star arrangement. This lends support to the picture that the ground state can be a VBC The hexagonal-VBC 0 , in contrast to the HVBC, lacks reflection symmetries about crystallographic axes; thus its symmetry group is reduced to C 6 . It has 12 classes of bonds. (c) The HVBC α has reduced (2π/3) rotation symmetry but preserves reflection symmetry; thus its symmetry group is C 3v . It has 14 classes of bonds.
with a 12-site unit cell capturing some modulation. In total, nine symmetry distinct VBCs with a 12-site unit cell can occur; see figures 3 and 4 for their NN patterns. In particular, the SVBC state (figure 3a) was argued in [52] to occur as an instability of the U(1) Dirac SL and to be consequently stabilized as the ground state of the NN spin-1/2 QHAF. Numerical studies using projected wave functions have shown this proposal to be incorrect and have also established the stability of the uniform RVB SL w.r.t. dimerizing into the SVBC state [32,34,35]. Furthermore, a recent QDM study [55] found the VBC 3 state (figure 4c) to be a competing ground state and a DMRG study [41] concluded that the DVBC state (figure 4b) is close by in a generalized parameter space. In section 3, we study the possibility of a ground state realization of VBC 3 and DVBC states numerically, within the framework of projected wave functions. In fact, we perform this study for all 12-site unit cell VBCs. ‡ The points P are not centers of inversion (π-rotation) symmetry, as has been mismarked in figure 1(a) of [55], which corresponds to the HVBC 0 state in the present work.  Figure 6: (a) The HVBC α 0 has only a reduced rotation symmetry (2π/3); thus in contrast to HVBC α its symmetry group is reduced to C 3 . It has 24 classes of bonds. (b) The 36-diamond-VBC has two perpendicular axes of reflection symmetry, thus giving rise to C 2v symmetry, with 19 classes of bonds. (c) The 36-VBC 3 has only π rotation symmetry; thus its symmetry group is C 2 . It has 36 classes of bonds.

36-site unit cell VBCs
The building blocks of the kagomé lattice can take nontrivial forms such as a 2 √ 3 × 2 √ 3 expansion of the elementary 3-site unit cell, thus giving rise to a tilted 36-site unit cell. It was shown in [43,44] that such a construction maximizes the density of hexagons on which dimer resonances occur, thereby lowering the energy. This 36-site unit cell VBC was studied numerically using series expansion [45,46] and MERA [47], which found it to be a good approximation to the ground state of the NN spin-1/2 QHAF. Similar conclusions were also obtained from a QDM study [49,50,55]. Motivated by these findings we classify all 36-site unit cell VBC patterns on the kagomé lattice, which leads to the identification of nine symmetry distinct VBCs; see figures 5, 6 and 7 for their NN patterns.
In our previous work [35] we studied the HVBC state (see figure 5a) by using projected wave functions and found it to be higher in energy compared to the gapless SLs. However, the symmetry of the VBC identified in QDM studies [49,50,55] is that of the HVBC 0 state (see figure 5b), which has a lower symmetry compared to the HVBC state. In section 3, we study the possibility of a ground state realization of the HVBC 0 state for the NN and NNN spin-1/2 QHAF.  Figure 7: (a) The 36-VBC 1 possesses only a single axis of reflection symmetry, which bisects the sides of the shaded hexagon; consequently, its symmetry group is C 1v . It has 38 classes of bonds. (b) The 36-VBC 1 has the same symmetry as 36-VBC 1 , but its reflection symmetry axis passes through a vertex of the shaded hexagon; we shall denote its symmetry group as C 1v , to distinguish it from that of 36-VBC 1 . It has 36 classes of bonds. Note that the 36-VBC 2 (which has no symmetry at all) has not been drawn. Its symmetry group is just the identity E. Consequently, it has 72 distinct classes of bonds.

General remarks on the VBC classification
It is worth mentioning that this VBC classification (for a given unit cell) is based on very general considerations of symmetry only and hence is not dependent on the formalism in which one studies these phases. In principle, it is possible to translate its construction from one language (e.g. QDM) into another (e.g. Schwinger fermions or bosons) for a VBC with a given symmetry. Moreover, within a given framework there can be different ways of constructing wave functions for a given VBC, consistent with its symmetry group. Firstly, one can add amplitudes beyond NN, consistent with the VBC symmetry group. Since we will study these phases within a slave particle approach, one can construct at the naive level simple mean-field wave functions or go much beyond mean-field and include the effects of full projection. At a next level, it is possible to improve the wave function by applying the Hamiltonian operator on it a given number of times and considering an optimized linear superposition of these wave functions with the original projected wave function. It is also worth noting that this hierarchical sorting of VBCs in each fixed symmetry sector also greatly eases the numerical search for a possible VBC stabilization as the ground state of the spin-1/2 QHAF.

Numerical Results
We study the energetics of SL and VBC phases for the spin-1/2 QHAF using Gutzwiller projected fermionic wave functions with the variational quantum Monte Carlo technique. Our variational calculations are performed on clusters with 432 (i.e. 3 × 12 × 12) or 576 (i.e. 36 × 4 × 4) sites and mixed periodic-antiperiodic boundary conditions which ensure non-degenerate mean-field wave functions at half filling. The large size of the cluster ensures that the spatial modulations induced in the observables by breaking of rotational symmetry (due to mixed boundary conditions) remain smaller than the uncertainty in the Monte Carlo simulations.
Among the class of NN fully symmetric and gapless SLs, the U(1) Dirac SL ([0, π] SL) has the lowest energy for the NN spin-1/2 QHAF. On a 432-site cluster its energy per site is E/J 1 = −0.428 63 (2) (1) for the uniform RVB SL [35]. Upon inclusion of NNN hopping amplitudes, one gets the extended U(1) Dirac SL or the extended uniform RVB SL, which are labeled by one additional flux through a plaquette of the type '234' in figure 1a, the flux through the other triangular plaquette formed by NNN bonds only is then fixed. Hence, the extended Dirac SL can be either the [0, π; π, 0] or the [0, π; 0, π] SL and analogously the extended uniform RVB SL can be either the [0, 0; π, π] or the [0, 0; 0, 0] SL [35]. For the NN spin-1/2 QHAF these extended SLs have a slightly lower energy, but they perform much better for the J 1 − J 2 spin-1/2 QHAF, see figure 9.

Results on the stability of gapless SLs towards VBC perturbations
We carried out an extensive numerical study of the local and global stability of the NN U(1) Dirac and uniform RVB SL toward dimerizing into all 6-, 12-and 36-site unit cell VBCs. In cases where we did find dimerization with NN bond amplitudes, we added second NN bond amplitudes to the SL and VBC ansatz (consistent with symmetries), since this led to a significant gain in energy. Our main focus was on the CVBC (figure 1b), DVBC (figure 4b), VBC 3 (figure 4c) and HVBC 0 (figure 5b) states, since these have been identified as ground states of the spin-1/2 QHAF in other studies. We perform our analysis by first fixing a background flux corresponding to the SL liquid whose stability we wish to study. Then, we introduce an amplitude modulation of χ ij consistent with the symmetries of the VBC, i.e. bonds belonging to the same class (color/line marking in figures 1b, 3, 4, 5, 6 and 7) have the same amplitude (χ λ ), which is set to different values for different classes. Starting from an arbitrary unbiased point (χ λ 's) in the variational space we perform an optimization of the wave function to obtain the lowest energy state [62,63].
3.1.1. The case of the U(1) Dirac SL For the NN spin-1/2 QHAF, the variation of parameters and energy in the Monte Carlo optimization for the four competing VBCs   (regarded as a dimerization of the U(1) Dirac SL) mentioned above is given in figure 8. As can be clearly seen, the energy converges neatly to the reference value of the U(1) Dirac SL, and all the parameters converge to χ λ = 1 (within error bars) after averaging over a sufficient number of converged Monte Carlo steps; thus the translation symmetry associated with the SL is restored. In fact, we performed these calculations for all 6-, 12-and 36-site VBCs and found that in each case the U(1) Dirac SL is stable towards opening a gap and destabilizing into any of these VBCs. This remarkable stability (for all VBCs) is also preserved upon addition of a NNN (J 2 ) super-exchange of hopping amplitudes and consequently a lower energy that is seen from the fact that the level crossing or the onset of VBC order is shifted from J 2 ≈ −0.09 [35] for HVBC to J 2,c ≈ −0.045 for HVBC 0 state. Thus our results still point to a gapless ground state for J 2 −0.045, which is along the lines of our previous work [35,36].

Conclusions and discussions
In this paper, we enumerated all 6-, 12-and 36-site unit cell VBCs based on symmetry considerations alone and subsequently investigated the possibility of stabilizing any of these VBCs in the NN and NNN spin-1/2 QHAF on a kagomé lattice. We found that the U(1) Dirac SL is remarkably robust toward dimerizing into any of these VBCs, for both the NN and NNN spin-1/2 QHAF. However, the uniform RVB SL dimerizes into a 36-site unit cell VBC, which becomes the lowest in energy on addition of a very weak FM coupling, J 2,c ≈ −0.045. Our systematic and thorough numerical investigation brings us to the conclusion that, at least within the Schwinger fermion approach to the spin model, the U(1) Dirac SL has the best variational energy for J 2 −0.045. The conflict between our results, which point to a gapless ground state in this region, and those obtained by exact diagonalizations and DMRG calculations, which instead suggested the presence of a fully gapped spectrum, remains open and deserves further investigation. One possible direction would be to include vison dynamics in the projected wave functions [64], which may be necessary to capture topological order faithfully. Another step would be to improve our variational wave functions based on the application of a few Lanczos steps [65] and then perform an approximate fixed-node projection technique. The possibility that an unconventional VBC breaking time-reversal symmetry is stabilized as the ground state cannot be ruled out [55]. Finally, we mention that VBC order might also set in via confinement transitions of the Z 2 SLs [66], this remains to be investigated numerically.