# Magnetic order and paramagnetic phases in the quantum -- honeycomb model

## Abstract

Recent work shows that a quantum spin liquid can arise in realistic fermionic models on a honeycomb lattice. We study the quantum spin-1/2 Heisenberg honeycomb model, considering couplings , , and up to third nearest neighbors. We use an unbiased pseudofermion functional renormalization group method to compute the magnetic susceptibility and determine the ordered and disordered states of the model. Aside from antiferromagnetic, collinear, and spiral order domains, we find a large paramagnetic region at intermediate coupling. For larger within this domain, we find a strong tendency to staggered dimer ordering, while the remaining paramagnetic regime for low shows only weak plaquet and staggered dimer response. We suggest this regime to be a promising region to look for quantum spin liquid states when charge fluctuations would be included.

###### pacs:

75.10.Kt, 75.10.JmThe search for a quantum spin liquid phase in nature ever since has been a complicated task not only from experiment, but also from theory (1); (2). This stems from the fact that the properties of a quantum spin liquid are very peculiar: quantum fluctuations have to form a many-body singlet state without long-range correlations of any kind of operator. Since the notion of frustration in quantum spin systems which is the main resource of fluctuations to accomplish a magnetically disordered quantum state at zero temperature, studies of a plethora of spin Hamiltonians on different lattices tell us that even most frustrated quantum systems tend to establish some sort of long-range correlations such as seen in a valence bond solid phase: while local spin operator correlations rapidly fall off there, long-range dimer-dimer correlations spoil the phenomenology of a quantum spin liquid.

The remarkable numerical studies by Meng et al. (3) are a fortunate exception. In their work, they report on the first unambiguous discovery of a genuine spin liquid phase from a generic microscopic model. They consider the Hubbard model on the honeycomb lattice by Monte Carlo methods and find a spin-gapped phase at which shows no long-range correlations of any kind, neither charge density wave, superconductivity or even spin solid-type correlations such as that of a valence-bond crystal formation. Moreover, the study finds a clear excitation gap and no symmetry breaking of any lattice symmetry or parity and time-reversal , which already excludes chiral spin liquids and algebraic spin liquids.

One path of further understanding the magnetic quantum phases on the honeycomb lattice is the development of effective descriptions for the Hubbard model itself such as gauge theory (4) and slave boson mean field theory methods (5); (6). Another direction, however, which we pursue in this Letter is to analyze the Gutzwiller-projected Hubbard model on the honeycomb lattice. Specifically, we consider Heisenberg spin couplings up to third nearest neighbors labeled as the -- model. The motivation for this is two-fold. First, the -- model projects onto the square root fraction of the Hilbert space of single site occupancy where only spin modes are present, which enables us to analyze the model through methods designed for this setup. Second, in reverse, a detailed understanding of the -- phase diagram will eventually help to identify which aspects of possible quantum phases may be explained through spin fluctuations only and which may necessitate the effect of charge fluctuations.

We consider the Hamiltonian

(1) |

where the sums extend over nearest, second nearest and third nearest neighbors, respectively (see Fig. 1). This model is expected to describe a class of magnetic materials with a honeycomb lattice, one example being (7). and are given in units of . The solution of the classical -- model has been known for a long time (8); (9). For small the system is Néel-ordered, which is commensurate with the honeycomb lattice and preserves the sublattice degrees rotational symmetry. For sufficiently low and beyond a threshold of , the system resides in a spiral phase. For high and , it is energetically preferred to order in a collinear phase where spins along zigzag chains align ferromagnetically while neighboring zigzag chains exhibit antiparallel spin orientation (there are three degenerate collinear configurations).

In order to obtain an adequate quantum phase diagram of the model defined in (1), there are not many suitable methods available, some of which predominantly focused on the - line. Exact diagonalization (ED) studies are very helpful, as arbitrary -point correlation functions can be computed in principle. However, except of the valence bond crystal domain where the dimer Hilbert space projection (10) is justified (11), ED cannot reach sufficient system sizes to adequately determine all phase regimes (12); (13); (11). Linear spin wave theory gives a qualitative tendency where the quantum corrections lead to, but is still strongly biased towards the classical limit (12); (13); (14). Similar shortcomings apply to Schwinger boson mean field theory (15); (13). A promising direction has recently been given by variational Monte Carlo (VMC) schemes, where energies of the antiferromagnetic (AFM) ordered trial state as well as several spin liquid (solid) candidates have been compared for , (16). There, a transition point of is found from the AFM ordered phase to the spin liquid phase which breaks no rotational invariance. It is followed by a transition to a dimerized spin solid phase at which breaks rotational invariance of the lattice. It is hence likely that the disordered phase itself, already for the - line and as such definitely for the full phase diagram with finite , contains different types of paramagnetic phases, which we will investigate in the following.

We employ the pseudo-fermion functional renormalization group (PFFRG) (17); (18) to identify magnetically ordered and paramagnetic phases and compute the groundstate magnetic susceptibility to determine the magnetic order parameter in the ordered phases as well as the dominant fluctuation profile in the paramagnetic phases. Our starting point is the pseudofermion representation of spin-1/2 operators , (, ) with the fermionic operators and and the Pauli-matrices . This representation enables us to apply Wick’s theorem leading to standard Feynman many-body techniques. Quantum spin models are inherently strongly coupled models, requiring an infinite self-consistent resummation theory. In this context FRG (19); (20); (21); (22); (23); (18); (24) provides a systematic summation in different interaction channels by generating equations for the evolution of all one-particle irreducible -particle vertex functions under the flow of an IR frequency cutoff . In order to reduce the infinite hierarchy of equations to a closed set, we restrict the computation to the full set of one-loop parquet diagrams and their vertex corrections up to infinite order, and additionally include certain two-loop contributions that are essential to induce a self-consistent resummation procedure exceeding the perturbative limit (25); (18). The parquet diagrams include graphs that favor magnetic order and those that favor disorder tendencies such that in total the method provides an adequate treatment of order and disorder fluctuations (26); (27); (18). From the two-particle scattering vertex we obtain the spin-spin correlation function (or spin susceptibility), which is the central outcome of our PFFRG. A magnetic ordering instability is initially signalled by a strong rise of the susceptibility associated with this order at some finite scale of . The onset of spontaneous long-range order is signalled by a sudden breakdown of the smooth flow (as shown for in Fig. 2). In the momentum-resolved magnetic susceptibility , where is defined along components in the extended Brillouin zone (BZ), the different magnetic ordering patterns manifest as peaks depicted in Fig. 2. Note that due to the two-atomic unit cell, the spin susceptibility does not have the periodicity of the first BZ but rather of the extended (second) BZ. Adding small dimer-field perturbations to the Hamiltonian, we are able to calculate dimer responses for a given dimer configuration (Fig. 2). For the present study our PFFRG algorithm internally deals with spin-spin correlations up to a length of lattice constants corresponding to a correlation area of sites which provides an adequate description of the model.

Fig. 1 shows the quantum phase diagram of the -- model (1). For dominant the system displays AFM order which persists longer against for finite as cooperates with . Increasing (for not too large ) we observe a melting of the order and the appearance of a rather large paramagnetic region. Above the system is characterized by presumably weak magnetic order and very small ordering instability scales which are hard to resolve numerically. However, as we enter this region by increasing , we observe the appearance of clean magnetic response peaks which we interpret as the onset of weak magnetic order. From the peak positions in -space we divide this region into a collinear ordered phase for large and a spiral ordered phase for smaller (Fig. 1). Throughout parameter space, the spiral phase is very close to -Néel order on both honeycomb sublattices except for a small region near the axis where the wave vector deviates from commensurability (Fig. 2). From the perspective of the degenerate classical spiral phase, this corresponds to a pinning of the order due to quantum fluctuations. A pronounced jump of the leading susceptibility wave-vector is seen as we cross the transition between AFM and C-AFM order, pointing to a first order transition. The observations are consistent with the classical result except for the fact that the transition between AFM and C-AFM order is shifted towards higher compared to due to quantum corrections included in our calculation.

We now focus on selected cuts through parameter space. We first investigate how the fluctuation profile changes along the line (Fig. 2). For small AFM order manifests itself in peaks at the corners of the extended BZ and a characteristic instability breakdown of the flow (lower line of Fig. 2). As we increase , the AFM peaks rapidly decrease and from the disappearance of unstable flow behavior we estimate the transition to be at . Inside the paramagnetic phase, such as at depicted in Fig. 2, no clear peak structure is visible and the susceptibility flow remains stable up to . Around spiral order peaks emerge at wave vectors slightly shifted from the commensurate positions towards larger . This feature is consistent with a large- expansion (14) which allows to select specific wave vectors out of a classical manifold of degenerate momenta. Upon further increasing , the peak positions approach commensurability, i.e., the sublattices effectively decouple and individually exhibit -Néel order. During the flow, these susceptibility peaks emerge at very small -scales, giving us indication of a very weak magnetization. While the transition between paramagnetism and spiral order is a bit smeared out at small , the onset of response peaks occurs more abruptly at larger .

To resolve more information about the correlations in the disordered phase, we compute the staggered and plaquet dimer response, i.e., the (dimensionless) factor of amplification of an external dimer-field perturbation exerted on the system. In doing so we can distinguish between parameter regimes of different dimer fluctuation strength. As we sweep through the paramagnetic phase at we find that the staggered dimer response is dominant for higher while plaquet and staggered dimer response compete for lower (Fig. 2). The absolute dimer response amplitudes are smaller for lower . This is in qualitative agreement with VMC (16) as well as with ED studies (11). It supports the view that if at all the system may form a spin liquid phase around this domain which is also the parameter regime related to the honeycomb Hubbard model from a strong coupling expansion. There, charge fluctuations which are neglected in (1) may be sufficient to destroy the comparably low dimer ordering tendency.

In addition, we investigate parameter lines varying for intermediate through the disordered regime. The fluctuation profiles and dimer responses for the paramagnetic regime at finite are shown in Fig. 3. We see that at small the fluctuation profiles for and are very similar but differ more with increasing until eventually the line leads into the AFM ordered phase while collinear order emerges on the line. The paramagnetic phase shows rather complicated susceptibility profiles which have lost most signature of ordering fluctuations. A typical feature is the ringlike shape as seen e.g. at . An intuitive reason for a quantum disordered phase is already indicated from the classical limit where the point is tricritical with three competing ordering tendencies. From the dimer responses along the and line we find that rather independent of , staggered dimer ordering tendency is more efficiently established by larger .

Note added. When this manuscript was completed we became aware of an independent work providing an analysis of model (1) through joint mean field and exact diagonalization techniques (28). Several similar findings show a good correspondence of both approaches.

###### Acknowledgements.

We thank S. Capponi, B. Clarke, A. Hackl, M. Hastings, A. Läuchli, C. Lhuillier, A. Muramatsu, K. P. Schmidt, R. Singh, S. Trebst, S. Wessel, and, in particular, P. Wölfle for discussions. RT thanks the KITP workshop ’Disentangling Quantum Many-body Systems: Computational and Conceptual Approaches’ for hospitality. JR is supported by DFG-FOR 960. RT is supported by DFG-SPP 1458/1 and the Humboldt Foundation.### References

- P. W. Anderson, Science 235, 1196 (1987).
- L. Balents, Nature 464, 199 (2010).
- Z. Y. Meng, T. C. Lang, S. Wessel, F. Assaad, and A. Muramatsu, Nature 464, 847 (2010).
- M. Hermele, Phys. Rev. B 76, 035125 (2007).
- F. Wang, Phys. Rev. B 82, 024419 (2010).
- A. Vaezi and X. G. Wen, arXiv:1010.5744.
- A. A. Tsirlin, O. Janson, and H. Rosner, Phys. Rev. B 82, 144416 (2010).
- E. Rastelli, A. Tassi, and L. Reatto, Physica B 97, 1 (1979).
- S. Katsura, T. Ide, and T. Morita, J. Stat. Phys. 42, 381 (1986).
- D. Poilblanc, M. Mambrini, and D. Schwandt, arXiv:0912.0724.
- H. Mosadeq, F. Shahbazi, and S. A. Jafari, arXiv:1007.0127.
- J. B. Fouet, P. Sindzingre, and C. Lhuillier, Eur. Phys. J. B 20, 241 (2001).
- D. C. Cabra, C. A. Lamas, and H. D. Rosales, arXiv:1003.3226.
- A. Mulder, R. Ganesh, L. Capriotti, and A. Paramekanti, Phys. Rev. B 81, 214419 (2010).
- A. Mattson, P. Fröjdh, and T. Einarsson, Phys. Rev. B 49, 3997 (1994).
- B. K. Clark, D. A. Abanin, and S. L. Sondhi, arXiv:1010.3011.
- J. Reuther and P. Wölfle, Phys. Rev. B 81, 144410 (2010).
- J. Reuther and R. Thomale, Phys. Rev. B 83, 024402 (2011).
- C. Wetterich, Phys. Lett. B 301, 90 (1993).
- T. R. Morris, Int. J. Mod. Phys. A 9, 2411 (1994).
- C. Honerkamp, M. Salmhofer, N. Furukawa, and T. M. Rice, Phys. Rev. B 63, 035109 (2001).
- R. Hedden, V. Meden, T. Pruschke, and K. Schönhammer, J. Phys.: Condes. Matter 16, 5279 (2004).
- M. Salmhofer, C. Honerkamp, W. Metzner, and O. Lauscher, Prog. Theor. Phys. 112, 943 (2004).
- J. Reuther, P. Wölfle, R. Darradi, W. Brenig, M. Arlego, and J. Richter, Phys. Rev. B 83, 064416 (2011).
- A. A. Katanin, Phys. Rev. B 70, 115109 (2004).
- P. W. Anderson, Phys. Rev. 86, 694 (1952).
- J. B. Marston and I. Affleck, Phys. Rev. B 39, 11538 (1989).
- A. F. Albuquerque, D. Schwandt, B. Hetenyi, S. Capponi, M. Mambrini, and A. M. Läuchli, arXiv:1102.5325.