Phase transitions of the ionic Hubbard model on the honeycomb lattice

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.

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.

Results
The strongly correlated honeycomb lattice with staggered potential. We consider the ionic Hubbard model on the honeycomb lattice (see inset in where c { is (c is ) creates (destroys) an electron with spin s at site i. t is the hopping amplitudes of fermions over nearest-neighbor sites, and we set t~1:0 as the unit energy. U (Uw0) is the amplitude of the onsite repulsive interaction, and D 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 m is fixed so that the average occupancy is (hn A izhn B i)=2~n~1{d, where d is the hole density.
We begin with the tight-binding Hamiltonian with staggered potential on the honeycomb lattice, corresponding to that the interaction U~0 in the ionic Hubbard model. After the fourier transformation, we can get the dispersion of the free electrons, where In this system, there are two bands, and the energy gap of the two bands DE~2D. From the tight binding model analysis above, we can learn that the system can be adjusted to various phases: a semimetal when staggered potential D~0 and hole doping density d~0, a band insulator when staggered potential D=0 and hole doping d~0, and a normal metal when staggered potential D=0 (or D~0) and hole doping d=0. In this paper, we mainly study the correlation effects in the band insulator and hole doping band insulator.
Phase diagram of the ionic Hubbard model. 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 U for half filling and hole doping at staggered potential D~0:4 and temperature T~1=20 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 U~0, the system is band insulator and normal metal at half filling and hole doping cases, respectively. With the increase of the interaction U, 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 U~U AF : 4:2. Below the critical interaction U AF , 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 DE ; vDE : . 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 U. At the hole density d~0:03, the phase transition occurs at critical interaction U~U1 e 4:1. 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 d~0:1, the system is paramagnetic metal below a critical interaction UvU1 e 3:8, ferromagnetic metal above another critical interaction UwU2 e 5:9, and antiferromagnetic metal between those two interaction U1vUvU2.
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 d~0:1, the phase transition of paramagnetic metal to antiferromagnetic metal is at U~U1 e 3:6 (U~3:8 for 6-site cluster), and the antiferromagnetic metal to ferromagnetic metal is at U~U2 e 5:8 (U~5:9 for 6-site cluster). We describe below in details of the spectral and magnetic properties that lead to this diagram.
Local quantities and spectral properties for the half filling case d~0. 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 d~0 for different values of the staggered potential D with the average occupancy n~1. 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 D~0:4 and d~0. Here we will give more detailed description on the results for different values of staggered potential D.
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 dn, double occupancy D occ , staggered magnetization M S , and uniform magnetization M F . 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 s hn ms i for a~A and B, N c is the site numbers of the cluster. We also calculate the double occupancy defined by The staggered magnetization and uniform magnetization are defined as and respectively, where the sublattice magnetization is calculated as m az~2 N c X m[a (hn m: i{hn m; i) for a~A and B. The results for the staggered charge density dn and the double occupancy D occ as a function of interaction U for temperature T~1=20 are shown in Figs. 2(a) and 2(b). Due to the staggered on-site potential, dn is always nonzero, even thought the Hubbard U tries to suppress it. dn decreases monotonically as a function of U, and shows no discontinuity at U AF . In the weak interaction region, the electrons prefer to gather on the lower potential sublattice B. The system experiences an imbalance between the two sublattice, resulting in higher double occupancy, compared with the Hubbard model when D~0 and a nonzero staggered charge density. Such tendencies become stronger as D grows. In the ionic limit D §t, it is energetically favorable that all the electrons are in the sublattice B, producing unity of the staggered charge. As U 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 M S and uniform magnetization M F as a function of interaction U for temperature T~1=20 respectively. For a given D, there exists a threshold value U AF at which the staggered magnetization turns on with a jump. Both the value of the U AF and the amplitude of the jump in M s are decreasing functions of D. In the half filling case d~0, the uniform magnetization M F is almost zero, independent of the staggered potential D and interaction strength U.
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 s is the spin, a~A,B, and v is measured from the chemical potential m. The density of states are derived from the imaginary time Green's function G(t) using maximum entropy method 43 . The local density of states are shown in Fig. 3 for several values of U at staggered potential D~0:5. For a quantitative analysis of the gap around a Fermi level, we investigate the spectral gap DE : and DE ; 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 U for D~0:5 at temperature T~1=20. In the noninteracting system U~0, the local density of states can be computed analytically and is composed of two bands which are separated by a band gap D 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 Local quantities and spectral properties for the hole doping case d=0. In this section, we now turn to study the phase transitions of the hole doping (d=0) band insulator as a function of the interaction on the honeycomb lattice. The average occupancy nv1, 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 D~0:4 and the temperature T~1=20, the phase diagram for different values of hole doping d is displayed in Fig. 1. For small hole doping d, the system goes a phase transition from paramagnetic metal to antiferromagnetic metal at a critical interaction U1. At finite hole doping d above a critical value, the system has three phases as the interaction U varying, paramagnetic metal when UvU1, antiferromagnetic metal when U1vUvU2, and ferromagnetic metal when UwU2. Now we discuss how we use the single particle spectral and other local qualities to determine the phase boundaries U1 and U2.
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 r s (v) for both of the spin components at hole doping d~0:1 and staggered potential d~0:4. The single particle density of states r s (v) 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 U~1:0, 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 U~4:5, 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 U~6:0, 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 A(k,v~0) very much. The k-resolved spectral weight can be defined as which are the maxima of the spectral weight at zero temperature as a function of k. In Fig. 6 we present A(k,v~0) for the three different phases at hole doping d~0:1 and staggered potential D~0:4. In the hole doping case, the Fermi surface A(k,v~0) are six rings in the K and K 0 points located at the corners of the hexagon. When the interaction is small U~1:0, corresponding to paramagnetic metal, the Fermi surface is only weakly renormalized compared to the case of the interaction U~0 (Figs. 6(a1)(a2)). In the intermediate interaction region U~4:5, corresponding to antiferromagnetic metal, the distribution of   quasi-particle spectral near the K and K 0 points became anisotropic in different directions in each corners (Figs. 6(b1)(b2)). In the large interaction region U~6:0, corresponding to ferromagnetic metal, the Fermi surface are strongly renormalized at each corners, the peak of A(k,v~0) are broadened (Figs. 6(c1)(c2)).
In addition to the spectral properties, we also calculate four local quantities: the staggered charge density dn, the double occupancy D occ , the staggered magnetization M s , and the uniform magnetization M F obtained from Eqs.(4), (5), (6) and (7), respectively. Figs. 7(a)(b) show the staggered charge density dn and the double occupancy D occ as a function of U for hole doping d~0:05, d~0:06 and d~0:1 at staggered potential D~0:4. 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 M S and the uniform magnetization M F as a function of interaction for hole doping d~0:02, d~0:06 and d~0:1 at staggered potential D~0:4. For small hole doping, such as d~0:02, the system shows a phase transition from paramagnetic metal to antiferromagnetic metal in which the staggered magnetization M S turns on a finite value at U~U1, and the uniform magnetization M F is zero all the time. When the hole doping above a critical value, such as d~0:06 and d~0:1, the magnetic properties shows dramatic changes at the phase boundaries U1 and U2. For small U, UvU1, the magnetic order is not favored. For intermediate U, U1vUvU2, there is a nonzero staggered magnetization M S . When the U is large, UwU2, the nonzero staggered magnetization M S is suppressed, the uniform magnetization M F becomes a nonezero value. Both of the staggered magnetization M S and the uniform magnetization M F increase with the increase of the hole doping d.

Discussion
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 U c , with the single particle gap decreases firstly, reaches a minimum at a critical interaction U c , then increases upturn, and the antiferromagnetic order gives a finite value above U c . 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 currents [44][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 2 , 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 simulated 29,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.

Methods
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