Cluster dynamical mean field theory of quantum phases on a honeycomb lattice
We have studied the ground state of the half-filled Hubbard model on a honeycomb lattice by performing the cluster dynamical mean field theory calculations with exact diagonalization on the cluster-impurity solver. Through using elaborate numerical analytic continuation, we identify the existence of a ‘spin liquid’ from the on-site interaction to (between and ) with a smooth crossover correspondingly from the charge fluctuation dominating phase into the charge correlation dominating phase. The semi-metallic state exits only at . We further find that the magnetic phase transition at from the ‘spin liquid’ to the Néel antiferromagnetic Mott insulating phase is a first-order quantum phase transition. We also show that the charge fluctuation plays a substantial role on keeping the ‘spin liquid’ phase against the emergence of a magnetic order.
pacs:71.10.-w, 71.27.+a, 71.30.+h, 71.10.Fd, 71.10.Pm
Quantum phase transition is a fascinating physics subject, which describes an abrupt change of the ground state of a quantum many-body system tuned by a non-thermal physical parameter, often accompanied with a novel quantum emergence phenomenon (1). As a canonical quantum phase transition, the Mott transition, from metallic to insulating state tuned by electronic Coulomb interaction, is one of the most celebrated and difficult problems in condensed matter physics (2). The resultant insulating state, namely Mott insulator, usually adopts spontaneous symmetry breaking in two and three spatial dimension to form a long-range antiferromagnetic (AFM) order to release the spin entropy due to localized electrons. Theoretically, the simplest model to capture such physics is the standard one-band half-filled Hubbard model. In the large Coulomb interaction limit this model reduces to a standard Heisenberg model with an AFM order in its ground state.
Nevertheless, an insulating ground state without any spontaneous symmetry breaking, namely spin liquid, may arise if there is frustration (3); (4). Actually the spin liquid is a genuine Mott insulating state in a sense that it is adiabatically separated from a band insulator. Spin liquid has been one of the most intriguing issues in condensed matter physics since it was introduced nearly forty years ago (3) and continuously in intense research since it was further proposed to be a parent phase to likely lead to high T superconductivity (5). However, a spin liquid had not been verified for a two-dimensional (2D) standard Hubbard or Heisenberg model until a recent quantum Monte Carlo simulation was done for correlated fermions on a honeycomb lattice (6). Through finite size extrapolation the simulation surprisingly shows that a spin liquid emerges between the semi-metallic and AFM Mott insulating states. This has stimulated a new surge of strong discussion on the nature of such a quantum phase and the related phase diagram (7); (8).
To help to clarify this issue, we performed the cluster dynamical mean field theory calculations for the half-filled Hubbard model defined on a honeycomb lattice as,
where is the electron creation (annihilation) operator with spin ( or ) at lattice site , represents the summation over the nearest neighbors, is the nearest neighbor hopping integral, and with the on-site Coulomb repulsion .
Our calculated results are schematically summarized in Fig. 1. By the calculations, we show that a disordered phase of ‘spin liquid’ exists from to , at which it transforms into the Néel AFM insulating phase via a first-order quantum phase transition.
The dynamical mean field theory (DMFT) maps a quantum lattice model onto a single lattice site, a quantum impurity, dynamically coupled to a self-consistently determined bath of free electrons that represents the rest of the lattice (9). Thus the DMFT fully considers local quantum dynamical fluctuations, much beyond conventional mean-field methods. The DMFT has substantially improved our understanding on the nonperturbative properties of correlated electron systems, particularly the Mott transition. The cluster dynamical mean field theory (CDMFT) is a natural extension of the DMFT to include the missed short-ranged spatial correlations through a proper replacement of a single-site impurity by a cluster of lattice sites, which is constructed to reflect the lattice symmetry and local lattice structure features (10); (11); (12). The CDMFT has been successfully applied to study a variety of ordered phases, and opens an avenue to directly study quantum phase transitions.
Here it should be addressed that for the (C)DMFT the thermodynamic limit is naturally taken from the outset through a self-consistent procedure (9); (12). As a nonperturbative approach to treat many-body correlation effects, the (C)DMFT works well in the whole coupling regime and becomes exact in the two contrary limits of both noninteracting and infinite-interacting cases. For finite-size quantum Monte Carlo simulations or exact diagonalizations, in contrast, the thermodynamic limit is extrapolated through finite-size scaling. Correspondingly, the CDMFT allows for spontaneous symmetry breaking, while the finite-size approaches have difficulty in finding a long-range order and related phase transition or underestimate ordered phases. Thus the CDMFT and finite-size approaches are complementary to each other, both of which together can give more conclusive results than they alone.
In CDMFT calculations, the target in the self-consistent procedure is to obtain a lattice imaginary frequency local Green’s function matrix (subscripts and being site indices of a chosen cluster) by assuming its self-energy matrix identified as the one of the corresponding cluster-impurity Green’s function matrix , derived from the Dyson equation. In order to study a ground state, we apply exact diagonalization rather than quantum Monte Carlo simulation to solve a cluster-impurity model (13); (14). Specifically, we employed the robust Krylov-Schur algorithm based SLEPc (15) to accomplish the large-scale sparse matrix diagonalization efficiently and stably (16). Here we particularly emphasize that to carry out an elaborate numerical analytic continuation from an imaginary frequency Green’s function onto a real frequency retarded Green’s function is crucial to unambiguously identify whether or not an energy gap exists at a small on-site interaction (17); (18); (19), by checking the density of states (DOS) equal to . In such a way, the energy gap resolution can be achieved as high as , which is hardly reached by other methods.
We first performed the CDMFT calculations for one-dimensional (1D) Hubbard model at half-filling with a cluster-impurity model respectively containing two and four lattice impurity sites, which serves as a benchmark for the further calculations. The calculated results are reported in Fig. 2, in quantitative comparison with the Bethe Ansatz exact solution (20). As we see, the two-site CDMFT result has been already in excellent agreement with the exact solution. Especially, by using the numerical analytic continuation, we can unambiguously identify a finite energy gap immediately develops once is nonzero. Thus the short-ranged spatial correlations play a dominant role in local spectral functions and properties even though the spatial correlations have long-ranged power-law behavior in 1D Hubbard model (21). It is also well-known that the quantum fluctuations are much stronger in one-dimension than in higher dimensions. Hence it is highly expected that the higher dimensional CDMFT results are more reliable and encouraging than the one-dimensional ones.
Figure 3(a) schematically shows the cluster-impurity model constructed for a honeycomb lattice. Such a model reflects sixfold rotational symmetry and an impurity site with a lattice coordination number of 3, which are the essential features of a honeycomb lattice. We also add a direct link between each pair of the nearest bath levels so that we can include the propagation of an electron from one impurity site through the outside of the cluster (the bath) to any other impurity site in the calculations.
When the on-site interaction , the half-filled Hubbard model on a honeycomb lattice reduces to a set of free Dirac fermions with a linear DOS around the Fermi energy. As shown in Fig. 3(b), the numerical analytic continuation can well repeat the exact DOS. Particularly, around the Fermi energy both are the strictly same even though the DOS is not differentiable at the Fermi energy.
In this study, we define the magnetization = on a site. Being bipartite, a honeycomb lattice can be divided into two sublattices and . If a Néel AFM state appears, will alternatively take positive and negative along with sublattices and , namely being staggered magnetization. As we see from Fig. 4, when , a paramagnetic solution of is stable, namely no spin polarization on each site, but over this solution is no longer stable and abruptly jumps over 0.2. On the other hand, when a staggered magnetization solution with is stable, namely a Néel AFM phase takes over, but below immediately plummets to zero. Between and , these two solutions coexist. Such a hysteresis behavior indicates this magnetic transition is a first-order quantum phase transition.
We calculated the DOS for a small on-site interaction with extreme caution through elaborate numerical analytic continuation (18). The calculated DOS are then plotted in Fig. 5. Similar to the case of 1D Hubbard model, what we find is that there is also a definite energy gap opening at the Fermi energy once is nonzero. In comparison with the case of by checking the enclosed area, it is further shown that the corresponding states nearby the Fermi energy are clearly moved away from the gap rather than pushed to the two sides of the gap. To be specific, for and , the energy gap is found as small as and , respectively. For a rather large , a relatively large energy gap opens with a substantial portion of states moved away from the Fermi energy into below and above , corresponding to the Hubbard band states. For a further larger than , the system transforms into the Néel AFM phase. The on-site spin degeneracy is then lifted. The spin-up (spin-down) resolved DOS at sublattice is the same as the spin-down (spin-up) resolved DOS at sublattice, as shown in Fig. 6.
Figure 7 shows the energy gap as a function of the on-site interaction , extracted from the calculated DOS. For , the energy gap increases very slowly with , and the function can be represented as . After , the energy gap increasing becomes fast with increasing. Between and , the -dependence of the energy gap shows a hysteresis behavior with a sudden change of , corresponding to the first-order quantum phase transition. Thus a nonzero definitely induces an energy gap and makes the system be in an insulating phase.
We also calculated the equal time spin-spin correlation functions and occupancy-occupancy correlation functions , with being a pair of the nearest sites, next nearest sites, and next next nearest sites, respectively. We find that all these correlation functions keep the sixfold rotational symmetry as the honeycomb lattice does before the magnetic transition. This excludes spin valence bond structures, spin anisotropic structures, and charge density waves. We can also exclude another possibility to break the translation symmetry that every six sites form a benzene-ring-like plaquette to further form a plaquette-singlet valence bond solid pattern in a honeycomb lattice. If such a pattern exists, the states nearby the Fermi energy will be pushed to the two sides of the gap rather than moved away. Thus this nonmagnetic insulating phase from to (between and ) can be classified as a ‘spin liquid’ phase in the sense that it is tuned on by the on-site interaction .
To understand the underlying physics, we examine the double occupancy, defined as on a site. The ground state energy per site of Hamiltonian (1) is a function of the on-site interaction . Its derivative is nothing but the double occupancy, namely . Thus the double occupancy directly describes a quantum phase transition tuned by . In Fig. 8, the -dependence of the likewise shows a hysteresis behavior between and . This means the energy level crossing in the ground state through the magnetic transition as increasing, being a characteristic of a first-order quantum phase transition (1).
The on-site interaction tunes or controls the Hamiltonian (1) through the double occupancy . Moreover, the localization degree of an electron, as well as the local correlation effect, can be quantitatively described by the double occupancy. At half-filling, is between 0.25 and 0, corresponding respectively to full delocalization and complete localization. In addition, the magnetic moment is obtained by .
It is commonly thought that the Mott transition is driven by the strong local correlation effect due to the on-site interaction , which is marked by a vanishing or very small double occupancy with a large local moment on a site (9); (22). In contrast, for the Hubbard model on a honeycomb lattice, a small can immediately induce a small energy gap opening to tune the system into an insulating phase with a large double occupancy namely large charge fluctuation, as shown in Fig. 8. The calculations show that the small- induced energy gap is a consequence of the interplay between the zero-DOS at the Fermi energy (Dirac Cone band) and local charge correlation, not a conventional correlation-driven Mott insulating gap. On the other hand, the correlation effect will become dominating nearby the magnetic transition. Actually Fig. 7 has shown that the energy gap becomes a linear function of after the transition, which is the canonical behavior of a correlation-driven Mott insulator. Thus the ‘spin liquid’ states nearby are of charge fluctuation dominating while those nearby are of charge correlation dominating, and corresponding to the ones found by the quantum Monte Carlo simulations reported in Ref. (6). Nevertheless our calculations show that a smooth crossover connects these two contrary parts. Meanwhile it is also indicated that the charge fluctuation plays a substantial role on keeping the ‘spin liquid’ phase against the emergence of an AFM order.
In summary, we have performed the cluster dynamical mean field theory calculations, allowing for the spontaneous symmetry breaking, to study the ground state of the half-filled Hubbard model on a honeycomb lattice. We find that a ‘spin liquid’ exists from to about , in which the system takes a smooth crossover correspondingly from the charge fluctuation dominating phase into the charge correlation dominating phase, then it further transforms into the Néel AFM Mott insulating phase via a first-order quantum phase transition.
Acknowledgements.We would like to thank Ning-Hua Tong for very helpful discussions. ZYL sincerely thanks the hospitality of International Center of Quantum Materials of Peking University, where this manuscript was finalized. This work is partially supported by National Natural Science Foundation of China and by National Program for Basic Research of MOST (2011CBA00112), China.
- Understanding Quantum Phase Transitions, edited by L. D. Carr (CRC Press, 2011); S. Sachdev, Quantum Phase Transitions (Cambridge University Press, 2000).
- N. F. Mott and R. Peierls, Proc. R. Soc. (London) Ser. A 49, 72 (1937); F. Gebhard, The Mott Metal-Insulator Transition: Models and Methods (Spring, 1997).
- P. W. Anderson, Mater. Res. Bull. 8, 153 (1973); P. Fazekas and P. W. Anderson, Philos. Mag. 30, 423 (1974).
- L. Balents, Nature 464, 199 (2010).
- P. W. Anderson, Science 235, 1196 (1987).
- Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Muramatsu, Nature 464, 847 (2010).
- F. Wang, Phys. Rev. B 82, 024419 (2010); Y. M. Lu and Y. Ran, Phys. Rev. B 84, 024420 (2011); B. K. Clark, D. A. Abanin, and S. L. Sondhi, Phys. Rev. Lett. 107, 087204 (2011); C. Xu, Phys. Rev. B 83, 024408 (2011).
- W. Wu, Y.-H. Chen, H.-Sh. Tao, N.-H. Tong, and W.-M. Liu, Phys. Rev. B 82, 245102 (2010); A. Liebsch, Phys. Rev. B 83, 035113 (2011).
- A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Rev. Mod. Phys. 68, 13 (1996); G. Kotliar and D. Vollhardt, Physics Today 57, 53 (2004).
- S. Moukouri and M. Jarrell, Phys. Rev. Lett. 87, 167010 (2001).
- G. Kotliar, S. Y. Savrasov, G. Palsson, and G. Biroli, Phys. Rev. Lett. 87, 186401 (2001).
- T. Maier, M. Jarrell, T. Pruschke, and M. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
- M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 (1994).
- C. A. Perroni, H. Ishida, and A. Liebsch, Phys. Rev. B 75, 045125 (2007).
- Vicente Hernandez, Jose E. Roman, and Vicente Vidal, SLEPc: A Scalable and Flexible Toolkit for the Solution of Eigenvalue Problems, ACM Transactions on Mathematical Software 31, 351 (2005); J. E. Roman, E. Romero, and A. Tomas, SLEPc home page: http://www.grycap.upv.es/slepc (2010).
- We rewrote the code from scratch in a fully parallel computimg structure to exploit the parallel framework of SLEPc and meanwhile be able to incorporate a variety of cluster-impurity model configurations.
- V. I. Anisimov, A. I. Poteryaev, M. A. Korotin, A. O. Anokhin, and G. Kotliar, J. Phys.: Condens. Matter 9, 7359 (1997).
- We adopted the Padé approximant algorithm for numerical analytic continuation (19). An effective temperature of is set for a small U, and for a large U.
- H. J. Vidberg and J. W. Serene, J. Low Temp. Phys. 29, 179 (1977).
- E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
- H. J. Schulz, Phys. Rev. Lett. 64, 2831 (1990).
- W. F. Brinkman and T. M. Rice, Phys. Rev. B 2, 4302 (1970).