RAPID COMMUNICATIONS PHYSICAL REVIEW B 89, 220408(R) (2014) Phase diagram of a frustrated spin- 1 2 J 1 - J 2 XXZ model on the honeycomb lattice P. H. Y. Li, 1, * R. F. Bishop, 1 , † and C. E. Campbell 2 1 School of Physics and Astronomy, Schuster Building, University of Manchester, Manchester M13 9PL, United Kingdom 2 School of Physics and Astronomy, University of Minnesota, 116 Church Street SE, Minneapolis, Minnesota 55455, USA (Received 8 April 2014; revised manuscript received 26 May 2014; published 24 June 2014) We study the zero-temperature (T = 0) ground-state (GS) properties of a frustrated spin-half J XXZ 1 –J XXZ 2 model on the honeycomb lattice with nearest-neighbor and next-nearest-neighbor interactions with exchange couplings J 1 > 0 and J 2 ≡ κJ 1 > 0, respectively, using the coupled cluster method. Both interactions are of the anisotropic XXZ type. We present the T = 0 GS phase diagram of the model in the ranges 0    1 of the spin-space anisotropy parameter and 0  κ  1 of the frustration parameter. A possible quantum spin-liquid region is identified. DOI: 10.1103/PhysRevB.89.220408 PACS number(s): 75.10.Jm, 75.10.Kt, 75.30.Kz, 75.30.Gw Frustrated spin-half (s = 1 2 ) antiferromagnets with nearest- neighbor (NN) J 1 > 0 and competing next-nearest-neighbor (NNN) J 2 > 0 exchange couplings on the honeycomb lattice have attracted a great deal of interest in recent years. These have included the two specific cases where both couplings have either an isotropic Heisenberg (XXX) form (see, e.g., Refs. [1–16], and references cited therein) or an isotropic XY (XX) form (see, e.g., Refs. [17–23]). Although the classical (s →∞) versions of these two models have identical zero-temperature (T = 0) ground-state (GS) phase diagrams [1,2], their s = 1 2 counterparts differ in significant ways. Furthermore, there is not yet a complete consensus on the GS phase orderings for either model in the range 0  κ  1 of the frustration parameter κ ≡ J 2 /J 1 . Whereas both classical (s →∞) models have N´ eel order- ing for κ<κ cl = 1 6 , the spin- 1 2 models both seem to retain N´ eel order out to larger values κ c 1 ≈ 0.2, consistent with the fact that quantum fluctuations generally favor collinear over noncollinear ordering. The degenerate family of spiral states that form the classical GS phase for all values κ>κ cl is very fragile against quantum fluctuations, and there is broad agreement that neither s = 1 2 model has a stable GS phase with spiral ordering for any value of κ in the range 0  κ  1. The most interesting, and also most uncertain, regime for both s = 1 2 models is when 0.2  κ  0.4. For the XXX model the N´ eel order that exists for κ<κ c 1 ≈ 0.2 is predicted by different methods to give way either to a GS phase with plaquette valence-bond crystalline (PVBC) order [6,7,10–14] or to a quantum spin-liquid (QSL) state [5,9,15,16] in the range κ c 1 <κ<κ c 2 ≈ 0.4. By contrast, for the XX model the N´ eel xy planar [N(p)] ordering that exists for κ<κ c 1 is predicted by different methods to yield either to a GS phase with N´ eel z-aligned [N(z)] order [19,23] or to a QSL state [17,20] in a corresponding range κ c 1 <κ<κ c 2 . There is broad agreement for both models that for (1 >) κ>κ c 2 there is a strong competition to form the GS phase between states with collinear N´ eel-II xy planar [N-II(p)] and staggered dimer valence-bond crystalline (SDVBC) forms of order, which lie very close in energy. The (threefold-degenerate) N´ eel-II states, which break the lattice rotational symmetry, are ones in which * peggyhyli@gmail.com † raymond.bishop@manchester.ac.uk NN pairs of spins are parallel along one of the three equivalent honeycomb directions and antiparallel along the other two. Some methods favor a further quantum critical point (QCP) at κ c 3 >κ c 2 , at which a transition occurs between GS phases with SDVBC ordering for κ c 2 <κ<κ c 3 , possibly mixed with N-II(p) ordering over all or part of the region, and N-II(p) ordering alone for κ>κ c 3 . The intriguing differences between the two models motivate us to consider the so-called J XXZ 1 –J XXZ 2 model that interpo- lates between them. It is shown schematically in Fig. 1(a) and is described by the Hamiltonian H = J 1  〈i,j 〉 ( s x i s x j + s y i s y j + s z i s z j ) + J 2  〈〈i,k〉〉 ( s x i s x k + s y i s y k + s z i s z k ) , (1) where 〈i,j 〉 and 〈〈i,k〉〉 indicate NN and NNN pairs of spins, respectively, and s i = (s x i ,s y i ,s z i ) is the spin operator on lattice site i . We shall study the T = 0 GS phase diagram for the spin- 1 2 Hamiltonian of Eq. (1) on the honeycomb lattice in the range 0    1 of the spin anisotropy parameter that spans from the XX model (with  = 0) to the XXX model (with  = 1), and in the range 0  κ  1 of the frustration parameter. Henceforth we put J 1 ≡ 1 to set the overall energy scale. We note that both exact diagonalization (ED) of small finite lattices and density-matrix renormalization group (DMRG) studies of the XX model in particular find it especially difficult to distinguish the N-II(p) and SDVBC phases in the regime κ>κ c 2 in the thermodynamic limit, N →∞, where N is the number of lattice sites. For this reason it is particularly suitable to use a size-extensive method such as the coupled cluster method (CCM) that works from the outset in the N →∞ limit. We first describe some key features of the CCM and refer the reader to Refs. [11,12,23–31] for more details. Any CCM calculation starts with the choice of a suitable model (or refer- ence) state |〉. Here we use each of the N(p), N(z), and N-II(p) states shown schematically in Figs. 1(b)–1(d). In order to treat each lattice site on an equal footing we passively rotate each spin in each model state, so that in its own local spin-coordinate frame it points downwards (i.e., along the local negative z axis). In these local spin coordinates every model state thus takes the universal form |〉 = | ↓↓↓ · · · ↓〉 and the Hamiltonian has to be rewritten accordingly. The exact GS energy eigenket, |〉, 1098-0121/2014/89(22)/220408(5) 220408-1 ©2014 American Physical Society