|Home | About | Journals | Submit | Contact Us | Français|
Many-body problem on the honeycomb lattice systems have been the subject of considerable experimental and theoretical interest. Here we investigate the phase transitions of the ionic Hubbard model on the honeycomb lattice with an alternate ionic potential for the half filling and hole doping cases by means of cellular dynamical mean field theory combining with continue time quantum Monte Carlo as an impurity solver. At half filling, as the increase of the interaction at a fixed ionic potential, we find the single particle gap decreases firstly, reaches a minimum at a critical interaction , then increases upturn. At , there is a band insulator to Mott insulator transition accompanying with the presence of the antiferromagnetic order. Away from half filing, the system shows three phases for the different values of hole density and interaction, paramagnetic metal, antiferromagnetic metal and ferromagnetic metal. Further, we present the staggered particle number, the double occupancy, the staggered magnetization, the uniform magnetization and the single particle spectral properties, which exhibit characteristic features for those phases.
The correlation effects in the honeycomb lattice systems have been extensively studied, which result in a number of exotic phenomena in both theory and experiment, such as the correlated electrons in the graphene1,2 and Silicene3,4,5, topological Mott insulator6 and quantum spin liquid7,8. Most of those studies are based on the standard Hubbard model, one of the most popular models in the strongly correlated system. For half-filling case, the electrons on the honeycomb lattice are described by a non-interacting massless Dirac fermion model with linear low energy dispersion relation. The system is semimetal, in which the Fermi surface are only six isolated points at the corners of the Brillouin zone. For the peculiar nature of the Fermi surface, the interaction effects can be suppressed by the low density of states in the Fermi level9,10,11,12. Away from the half-filling, the different behavior will arise in this system13. For example, at the or filling, the system shows many weak coupling instabilities to various ordered states, including spin density waves14, Pomeranchuk metal15, and p/d-wave superconductors16,17,18.
Recently, a new class of two dimensional materials has been found19,20,21, which is formed on a single layer honeycomb lattice consisting of alternating “M” and “N” orbitals with a level offset. Experimental results show that supports a unconventional superconductor19,22. The origin of superconductivity can be revealed based on the ionic-Hubbard model on the honeycomb lattice with the staggered lattice potential. The ionic Hubbard model, an extended version of the Hubbard model, was proposed to explain the neutral-ionic transition in the quasi-one-dimensional charge-transfer organic materials23. It has also been proposed to investigate the band insulator to Mott insulator transition, such as the one dimension system24,25 and two dimension square lattice system26. However, the charge dynamics with spins and the phase diagram of this model on the honeycomb lattice have not been studied. Moreover, the ionic Hubbard model on the honeycomb lattice can also be realized by cold atoms loaded in the optical lattices27,28,29,30,31,32, in which the on-site interaction, hopping amplitude, doping, and temperature can be fully controlled using Feshbach resonances, changing the lattice depth, changing the number of fermions, and varying the cooling time.
The dynamical mean field theory (DMFT)33 and its cluster extensions34,35 are powerful method to investigate the strongly correlated system, due to the efficient description of the quantum fluctuations. The cellular dynamical mean field theory (CDMFT)35 is one of the cluster extensions of DMFT, and the cluster is constructed in real lattice space. In contrast to a single site is chosen to construct the self-consistent equation in DMFT, the CDMFT picks up a cluster. This makes it is possible to include short range spatial fluctuations inside the cluster, which are important in the low dimensional systems. This method have been used to study the correlation effects on the honeycomb lattice and square lattice, such as Mott transition36,37, topological phase transition38 and charge order insulator transition39,40.
In this paper, we study the phase transitions of the ionic Hubbard model on the honeycomb lattice as a function of the hole doping and temperature. We adopt the CDMFT combined with the continuous time quantum Monte Carlo method (CTQMC)41,42. In order to determine the phase diagram, we calculate the staggered particle number, the double occupancy, the staggered magnetization, the uniform magnetization and the single particle spectral properties. At half filling, the system goes from a paramagnetic band insulator phase to an antiferromagnetic Mott insulator phase with the increase of the interaction. At small hole doping, the system has two phases, a paramagnetic metal for weak interaction and an antiferromangetic metal for large interaction. For finite hole doping above a critical value, the system shows three phases, a paramagnetic metal at weak interaction region, a antiferromagnetic metal at intermediate interaction region, then a ferromagnetic metal at strong interaction region.
We consider the ionic Hubbard model on the honeycomb lattice (see inset in Fig. 1). The system is composed of two alternating sublattices and . The Hamiltonian can be written as
where creates (destroys) an electron with spin at site . is the hopping amplitudes of fermions over nearest-neighbor sites, and we set as the unit energy. is the amplitude of the on-site repulsive interaction, and is a staggered one-body potential on the two sublattices in each unit cell, which is also called the “ionic" potential. The last term, the chemical potential is fixed so that the average occupancy is , where is the hole density.
We begin with the tight-binding Hamiltonian with staggered potential on the honeycomb lattice, corresponding to that the interaction in the ionic Hubbard model. After the fourier transformation, we can get the dispersion of the free electrons,
In this system, there are two bands, and the energy gap of the two bands . From the tight binding model analysis above, we can learn that the system can be adjusted to various phases: a semimetal when staggered potential and hole doping density , a band insulator when staggered potential and hole doping , and a normal metal when staggered potential (or ) and hole doping . In this paper, we mainly study the correlation effects in the band insulator and hole doping band insulator.
In this section, we summarize our main results of the ionic Hubbard model on the honeycomb lattice, deferring the details of how they were obtained to the following sections. The phase diagram as a function of interaction for half filling and hole doping at staggered potential and temperature obtained from the analysis using 6-site cluster is shown in Fig. 1. The results obtained using 8-site cluster are also shown to quantitatively see the cluster-size dependence. In the noninteracting limit , the system is band insulator and normal metal at half filling and hole doping cases, respectively. With the increase of the interaction , the system shows two phases for the half filling case, corresponding to the band insulator and the antiferromagnetic Mott insulator, and the two phases separate at the critical interaction . Below the critical interaction , the energy gap in the band insulator are the same for both spin components and decrease as the interaction increasing. In the Mott insulator phase, the single particle energy gap are different for both spin components, such as . And the Mott gap increase monotonously with the increase of the interaction.
For the small hole doping case, the system goes a phase transition from paramagnetic metal to antiferromagnetic metal when changing the interaction . At the hole density , the phase transition occurs at critical interaction . For finite hole doping above a critical value, there are three phases at different interaction, corresponding to paramagnetic metal, antiferromagnetic metal, and ferromagnetic metal. For example, at hole doping , the system is paramagnetic metal below a critical interaction , ferromagnetic metal above another critical interaction , and antiferromagnetic metal between those two interaction .
In Fig. 1, we also present the results for the 8-site cluster. In this case, the properties of this system are qualitatively same, but the phase boundary shifts a little, such as, in , the phase transition of paramagnetic metal to antiferromagnetic metal is at ( for 6-site cluster), and the antiferromagnetic metal to ferromagnetic metal is at ( for 6-site cluster). We describe below in details of the spectral and magnetic properties that lead to this diagram.
In this section, we firstly try to understand the correlation effects on the band insulator on the honeycomb lattice. We concentrate on the half-filling case for different values of the staggered potential with the average occupancy . In the noninteracting limit, the system prefers a band insulator phase, in which most of the electrons stay on a sublattice with lower potential, resulting in zero density of states in the Fermi surface. When the local interaction is turned on, the band insulator competes with the Mott insulator with one electron per lattice site. In Fig. 1, the phase diagram give us the results of the band insulator to Mott insulator at and . Here we will give more detailed description on the results for different values of staggered potential .
In order to examine how the system evolves with the variation of the local interaction, we firstly calculate the four local quantities: staggered charge density , double occupancy , staggered magnetization , and uniform magnetization . The staggered charge density and double occupancy are related with the charge fluctuation, the staggered magnetization and the uniform magnetization give us the information about the spin fluctuation. The staggered charge density is defined by the difference between the particle number densities at two sublattices,
where the sublattice number densities can be calculated as for and , is the site numbers of the cluster. We also calculate the double occupancy defined by
The staggered magnetization and uniform magnetization are defined as
respectively, where the sublattice magnetization is calculated as for and .
The results for the staggered charge density and the double occupancy as a function of interaction for temperature are shown in Figs. 2(a) and 2(b). Due to the staggered on-site potential, is always nonzero, even thought the Hubbard tries to suppress it. decreases monotonically as a function of , and shows no discontinuity at . In the weak interaction region, the electrons prefer to gather on the lower potential sublattice . The system experiences an imbalance between the two sublattice, resulting in higher double occupancy, compared with the Hubbard model when and a nonzero staggered charge density. Such tendencies become stronger as grows. In the ionic limit , it is energetically favorable that all the electrons are in the sublattice , producing unity of the staggered charge. As increasing, the energy cost of two electrons to stay in the same site becomes large, both the double occupancy and the staggered charge density decrease monotonically with the imbalance between the two sublattices become weaker. In the strong coupling limit, the staggered charge density is close to 0.
In Figs. 2(c) and 2(d) we plot the staggered magnetization and uniform magnetization as a function of interaction for temperature respectively. For a given , there exists a threshold value at which the staggered magnetization turns on with a jump. Both the value of the and the amplitude of the jump in are decreasing functions of . In the half filling case , the uniform magnetization is almost zero, independent of the staggered potential and interaction strength .
The local density of states provide more detailed information on the single particle properties. The spin-resolved single particle density of states are computed as follows
where is the spin, , and is measured from the chemical potential . The density of states are derived from the imaginary time Green’s function using maximum entropy method43. The local density of states are shown in Fig. 3 for several values of at staggered potential . For a quantitative analysis of the gap around a Fermi level, we investigate the spectral gap and for both spin components which are defined as the energy difference between the highest filled and lowest empty levels in the local density of states. Fig. 4 shows the spectral gap as a function of the interaction strength for at temperature . In the noninteracting system , the local density of states can be computed analytically and is composed of two bands which are separated by a band gap due to the staggered potential. For weak interaction, the local density of states are the same for both spin components. However the band gap around a Fermi level decreases monotonically with the increase of interaction in this region. For example, the density of states for two spin components at interaction and are shown in Figs. 3(a1)(a2) and Figs. 3(b1)(b2) respectively, the band gap at is smaller than . The local density of states displays a minimum spectral gap at a critical value of interaction , for up-spin component and for down-spin component. When above the critical value, the gap between the two bands in turn is enlarged, and the system goes to the Mott insulator phase. In this region, the local density of states for up-spin component and down-spin component show different behaviors (Figs. 3(c1)(c2). The Mott gap of the up-spin component is bigger than the down-spin component.
In this section, we now turn to study the phase transitions of the hole doping () band insulator as a function of the interaction on the honeycomb lattice. The average occupancy , resulting in finite density of states in the Fermi level. The system is a metal, and the low energy physics can be described by the Fermi liquid theory. When the interaction is turned on, there are many instabilities for the Fermi liquids, which is an enduring theme research in condensed matter physics. When the staggered potential and the temperature , the phase diagram for different values of hole doping is displayed in Fig. 1. For small hole doping , the system goes a phase transition from paramagnetic metal to antiferromagnetic metal at a critical interaction . At finite hole doping above a critical value, the system has three phases as the interaction varying, paramagnetic metal when , antiferromagnetic metal when , and ferromagnetic metal when . Now we discuss how we use the single particle spectral and other local qualities to determine the phase boundaries and .
Firstly, let us study the single particle density of states as a function of interaction for various hole doping. Fig. 5 shows the single particle density of states for both of the spin components at hole doping and staggered potential $. The single particle density of states are obtained from Eq.(8). With the increase of the interaction, the spectral weight is continuously transferred to the higher energy states, the chemical potential lies inside the lower band for both spin components all the time. When , in the paramagnetic metal region, the density of states for both spins are the same, and there are two the spectral peaks above and below the Fermi level (Figs. 5(a1) and 5(a2)). When , corresponding to the antiferromagnetic metal, the antiferromagnetic order sets in, making the density of states and gaps a little different for the two spin components (Figs. 5(b1) and 5(b2)). When , in the ferromangnetic metal region, the density of states for the two spin components are renormalized much and very different (Figs. 5(c1) and 5(c2)). In both the antiferromagnetic metal and ferromagnetic metal, one of the spectral peaks above the Fermi level is suppressed.
Besides the changes of the local density of states, the interaction will influence the momentum-resolved spectral density in the Fermi level very much. The -resolved spectral weight can be defined as
which are the maxima of the spectral weight at zero temperature as a function of . In Fig. 6 we present for the three different phases at hole doping and staggered potential . In the hole doping case, the Fermi surface are six rings in the and points located at the corners of the hexagon. When the interaction is small , corresponding to paramagnetic metal, the Fermi surface is only weakly renormalized compared to the case of the interaction (Figs. 6(a1)(a2)). In the intermediate interaction region , corresponding to antiferromagnetic metal, the distribution of quasi-particle spectral near the and points became anisotropic in different directions in each corners (Figs. 6(b1)(b2)). In the large interaction region , corresponding to ferromagnetic metal, the Fermi surface are strongly renormalized at each corners, the peak of are broadened (Figs. 6(c1)(c2)).
In addition to the spectral properties, we also calculate four local quantities: the staggered charge density , the double occupancy , the staggered magnetization , and the uniform magnetization obtained from Eqs.(4), (5), (6) and (7), respectively. Figs. 7(a)(b) show the staggered charge density and the double occupancy as a function of for hole doping , and at staggered potential . Although the system is metallic all the times, with the increase of the interaction, both of the two quantities decrease monotonically. Figs. 7(c)(d) show the staggered magnetization and the uniform magnetization as a function of interaction for hole doping , and at staggered potential . For small hole doping, such as , the system shows a phase transition from paramagnetic metal to antiferromagnetic metal in which the staggered magnetization turns on a finite value at , and the uniform magnetization is zero all the time. When the hole doping above a critical value, such as and , the magnetic properties shows dramatic changes at the phase boundaries and . For small , , the magnetic order is not favored. For intermediate , , there is a nonzero staggered magnetization . When the is large, , the nonzero staggered magnetization is suppressed, the uniform magnetization becomes a nonezero value. Both of the staggered magnetization and the uniform magnetization increase with the increase of the hole doping .
In this work, we have investigated the effect of on-site interaction and staggered ionic potential in a band insulator and doped band insulator on the honeycomb lattice based on the ionic Hubbard model. By means of the cellular dynamical mean field theory combing with continue time Monte Carlo method, we construct a phase diagram as a function of interaction and hole doping. At half filling, although the single particle spectral functions always posses a energy gap, the system shows a band insulator to Mott insulator transition at a critical interaction , with the single particle gap decreases firstly, reaches a minimum at a critical interaction , then increases upturn, and the antiferromagnetic order gives a finite value above . Away from half filing, many metallic phases with magnetic order are found, in order to exhibit characteristic features of the phases, the behavior of the staggered particle number, the double occupancy, the staggered magnetization, the uniform magnetization and the single particle spectral properties have bend studied. At small hole doping, the system goes a phase transition from a paramagnetic metal to an antiferromangetic metal with the increase of the interaction. For finite hole doping above a critical value, the system shows three phases, a paramagnetic metal at weak interaction region, a antiferromagnetic metal at intermediate interaction region, then a ferromagnetic metal at strong interaction region.
We get itinerant metals with spin density wave state which are an interesting class of materials where electrons show spin polarization or staggered spin polarization behavior. They have applications in spintronics as they can generate spin-polarized currents44,45,46. And the materials with the honeycomb lattice structure are very common, such as single layer graphene, silicene considered as the silicon-based counterpart of graphene, and monolayer molybdenum disulfide (ML-MDS), MoS, which play a vital role in nanoelectronics and nanospintronics. We hope that our study will motivate a research on along those lines and open up new possibilities in the area of spintronics. Moreover, with the development of the cold atom experiment, the honeycomb lattice have been simulated29,32,47,48, which can give us a platform to simulate and detect the phase transitions by loading ultracold atoms on the honeycomb optical lattices.
In order to study the ionic Hubbard model in honeycomb lattice which describes the correlation effects on the band insulator and the hole doped band insulator, the Cellular dynamical mean field theory are employed. The Cellular dynamical mean field theory is an extension of dynamical mean field theory, which is able to partially cure dynamical mean field theory’s spatial limitations. We replace the single site impurity by a cluster of impurities embedded in a self-consistent bath. The cluster-impurity problem embedded in a bath of free fermions can be written in a quadratic form,
where i and j are the coordinates inside the cluster-impurity, and the is the Weiss field. The effective medium is computed via the Dyson equation,
Within cellular dynamical mean field theory, the interacting lattice Green’s function in the cluster site basis is given by,
where are Matsubara frequencies, is the chemical potential and is the Fourier-transformed hopping matrix for the super lattice. In our analysis, the 6- and 8-site clusters in the inset of Fig. 1 are used to set up the cluster Hamiltonian. For the 6-site cluster case, the hopping matrix of the cluster can be written as follows ( the reduced Brillouin-zone),
where are the nearest-neighbor super-lattice vectors, is the lattice constant. In each iteration, in order to solve the effective cluster model and to calculate , we use the weak coupling interaction expansion continuous time quantum Monte Carlo method.
The CDMFT iteration procedure is summary as follows. Given a cluster self-energy , we can compute via Eq.(5)(6), then solve the effective cluster model and to calculate a new . Then us Eq.(5) again, we can get a new cluster self-energy . Repeat the procedure until the results are convergence.
The weak coupling interaction expansion continuous time quantum Monte Carlo method is efficient method to treat the impurity model. The method employs same tricks, which used to derived Feynman perturbation theory, to stochastically generate the partition function . In the interaction picture, , where is the time-ordering operator. The expansion of the partition function in power of reads
where and is the number of the sites of the cluster. The observable expectation value can be sampled during the Monte Carlo update. For example, the Green’s function
H. F. Lin acknowledges very helpful discussions with Y. H. Chen, F. J. Huang and Q. H. Chen. This work was supported by the NKBRSFC under Grants No. 2011CB921502, No. 2012CB821305, NSFC under Grants No. 61227902, No. 61378017, and No. 11434015, SKLQOQOD under Grant No. KF201403, SPRPCAS under Grant No. XDB01020300. We acknowledge the supercomputing center of CAS for kindly allocating computational resources.
Author Contributions H. F. L. performed calculations. H. F. L., H. D. L., H. S. T., W. M. L. analyzed numerical results. H. F. L., W. M. L. contributed in completing the paper.