Intermediate Phase in Interacting Dirac Fermions with Staggered Potential
Jingyao Wang, Lufeng Zhang, Runyu Ma, Qiaoni Chen, Ying Liang and, Tianxing Ma

TL;DR
This study uses quantum Monte Carlo simulations to uncover an intermediate metallic phase driven by electronic correlations in interacting Dirac fermions with staggered potential, revealing a complex phase diagram with potential experimental implications.
Contribution
The paper identifies a novel intermediate metallic phase in a model of interacting Dirac fermions, driven by correlations, which was not previously characterized.
Findings
Intermediate phase is metallic and driven by correlations.
Mott insulator phase exhibits antiferromagnetic order.
Phase diagram shows the intermediate state is robust and experimentally detectable.
Abstract
By performing exact quantum Monte Carlo simulations of a model of interacting Dirac Fermions with staggered potential, we reveal a novel intermediate phase where the electronic correlations drive a band insulator metallic, and at a larger interaction, drive the metal to Mott insulator. We also show that the Mott insulating phase is antiferromagnetic. A complete phase diagram is achieved by studying the phase transitions at large staggered potential and interaction strengths, which shows that the intermediate state is robust and occupies a large part of the phase diagram and that it should be more feasible to be detected experimentally.
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Intermediate phase in interacting Dirac fermions with staggered potential
Jingyao Wang
Department of Physics, Beijing Normal University, Beijing 100875, China
Department of Physics and Astronomy, University of Southern California, Los Angeles, CA 90089-0484, USA
Lufeng Zhang
School of Science, Beijing University of Posts and Telecommunications, Beijing 100876, China
Department of Physics, Beijing Normal University, Beijing 100875, China
Runyu Ma
Department of Physics, Beijing Normal University, Beijing 100875, China
Qiaoni Chen
Department of Physics, Beijing Normal University, Beijing 100875, China
Ying Liang
Department of Physics, Beijing Normal University, Beijing 100875, China
Tianxing Ma
Department of Physics, Beijing Normal University, Beijing 100875, China
Beijing Computational Science Research Center, Beijing 100193, China
Abstract
By performing exact quantum Monte Carlo simulations of a model of interacting Dirac Fermions with staggered potential, we reveal a novel intermediate phase where the electronic correlations drive a band insulator metallic, and at a larger interaction, drive the metal to Mott insulator. We also show that the Mott insulating phase is antiferromagnetic. A complete phase diagram is achieved by studying the phase transitions at large staggered potential and interaction strengths, which shows that the intermediate state is robust and occupies a large part of the phase diagram and that it should be more feasible to be detected experimentally.
I Introduction
Since the discovery of grapheneNovoselov et al. (2005); *Zhang2005; Castro Neto et al. (2009) and topological insulatorsHasan and Kane (2010); *RevModPhys.83.1057, Dirac fermions described by a honeycomb lattice, has enriched our knowledge of physics beyond Landau’s symmetry breaking theoryPeres (2010); Mourik et al. (2012); Grover et al. (2014). Landau’s theory of the Fermi liquid describes the interacting electrons of a typical metal as an ideal gas of noninteracting quasiparticles. This description is expected to fail for Dirac fermions due to their linearly dispersing bands and minimally screened Coulomb interactionsCastro Neto et al. (2009). Half-filled graphene hosts a Dirac fluid governed by relativistic hydrodynamicsCrossno et al. (2016); *Gallagher158. Inspired by recent experiments on twisted bilayer grapheneCao et al. (2018a); *Cao2018B; Yankowitz et al. (2019); Huang et al. (2019); *PhysRevB.101.155413, a fast-growing frontier of research has focused on the novel physics induced by correlation effects in interacting Dirac fermions.
Correlation effects play an essential role in many intriguing physical properties of modern science, touching upon topics ranging from unconventional superconductivityDagotto (1994); Lee et al. (2006), fractional quantum Hall effectXu et al. (2009); Bolotin et al. (2009), quantum spin liquidSorella et al. (2012); *PhysRevX.6.011029; Zhou et al. (2017); Liu et al. (2018), to metal-insulator phase transitionsMa et al. (2018); Paris et al. (2007). Those phenomena are all relevant to the ionic Hubbard modelHubbard and Torrance (1981), which contains the on-site Coulomb interactions and staggered potentials on bipartite lattices. Generally speaking, in a bipartite lattice like the honeycomb lattice, a broken inversion symmetry caused by an energy offset between the two sublattices, leads to a trivial band insulator at half filling, and the Coulomb interaction slows down the electrons or even localizes them in a Mott insulating phase, characterized by a spectral gap that opensJiang et al. (2019). Studies on the competition between the on-site Coulomb interactions and staggered potentials have witnessed extraordinary growth about possible exotic intermediate states between two or more competing phases.
This issue was very actively debated over more than a decadeEgami et al. (1993); Jiang et al. (2019); Hubbard and Torrance (1981); Resta and Sorella (1995); Batista and Aligia (2004); Niyaz et al. (1991); Rousseau et al. (2006); Kancharla and Dagotto (2007); Ebrahimkhas et al. (2015); Bag et al. (2015) before finally being settledMesser et al. (2015); Loida et al. (2017). The seminal work of dynamical mean-field theory (DMFT) studies on a square latticeGarg et al. (2006) suggests that an interaction-induced metallic phase exists in the intermediate coupling region. Subsequently, cellular DMFT simulations have found a bond order phaseKancharla and Dagotto (2007) in this region, while determinant quantum Monte Carlo (DQMC) calculations of conductivity indicate a metallic phaseParis et al. (2007); Bouadim et al. (2007). In addition to metallic and bond order insulating phases, various other phases depending on the lattice geometry have been proposed on other bipartite lattices, such as charge-density-wave insulatorKancharla and Dagotto (2007), superfluidNiyaz et al. (1991); Rousseau et al. (2006), semimetalEbrahimkhas et al. (2015) and half-metalBag et al. (2015). Recently, exciting progress on ultracold atom experiments has been made, and the ionic Hubbard model was realized in an optical honeycomb latticeLoida et al. (2017); Messer et al. (2015). Unfortunately, only the band insulating phase and Mott insulating phase were observed. Therefore, determining the existence of an intermediate phase or the nature of the intermediate phase is a subtle and largely open problem. Stimulated by the controversy and to motivate further experiments, in this paper, we focus on the ionic Hubbard model on a honeycomb lattice, which is a minimal model that includes both interactions and staggered potentials in a two-dimensional (2D) Dirac system. This model can not only be implemented in cold-atom systems but also can be realized on hydrogen grapheneBalog et al. (2010); moreover, the new family of 2D layered nitride materials Lix**MNCl (M=Hf, Zr) may be another candidateYamanaka et al. (1998); *PhysRevLett.94.217002.
Our simulations were completed by the DQMC method on a half-filled case. By varying the on-site interaction , staggered potential , lattice size and temperature, the bulk conductivity is calculated to determine either a metallic or an insulating phase, and finite-size scaling is implemented to detect the long-range antiferromagnetic (AFM) order in the thermodynamic limit. Our results reveal that an intermediate metallic state exists between the band and Mott insulators, and the AFM order appears after the second metal-insulator transition with increasing . The exact numerical method that we are using successfully captures all the phase transitions at large-enough staggered potential and interaction strength, and a complete phase diagram is achieved and shown in Fig. 1, which has several key differences relative to previous models. First, the intermediate state that we found is more robust and occupies a larger part of the phase diagram. For example, at small , the intermediate phase disappears quickly as increases in the ionic Hubbard model on a square latticeParis et al. (2007), while here, the intermediate state is robust up to =3.9. Second, the critical is in a reasonable range for experimental detection. The intermediate insulator state in the square lattice vanishes around =11Kancharla and Dagotto (2007) while it continues in the Haldane-Hubbard ModelVanhala et al. (2016). Furthermore, beyond previous results, we show a complete phase diagram where for large , the intermediate state disappears, and the system transitions from band insulator to Mott insulator directly.
II Model and method
The Hamiltonian for interacting Dirac Fermions with staggered potential is
[TABLE]
where , and represent the nearest-neighbor electron hopping amplitude, on-site Coulomb repulsion, and chemical potential. The electron density of the system is characterized by the chemical potential . is the operator that creates (annihilates) an electron with spin at site , and =. Specifically, is the staggered one-body potential between sites in A and B sublattices with opposite signs. It is also called the ionic potential. The band gap, , has a nonzero value as a result of breaking the symmetry between sublattice A and B. We set =1 as the default energy scale.
We adopt the exact DQMC methodBlankenbecler et al. (1981); *PhysRevB.40.506 to study the phase transitions in the model defined by Eq.(1). DQMC is a powerful and reliable tool to investigate strongly correlated electron systems. In the DQMC method, the partition function = is expressed as a path integral over the discretized inverse temperature over a set of random auxiliary fields. Then, the integration is accomplished by Monte Carlo techniques. The on-site interaction is decoupled by a Hubbard-Stratonovich (HS) transformation, which leads to a sum over the discrete HS field and leaves the Hamiltonian in a quadratic form in the fermion operators. The resulting quadratic form can be integrated analytically and becomes the Boltzmann weight, expressed as the product of the determinants of two matrices that depend on the HS spin variables. In the half-filled ionic Hubbard model on the honeycomb lattice, the system is sign-problem free on account of the particle-hole symmetry, which allows us to achieve large simulations to converge to the ground state.
With the aim of exploring the phase transitions of the system, we computed the -dependent DC conductivity, which is calculated from the wave vector q and the imaginary time -dependent current-current correlation functionTrivedi et al. (1996); Trivedi and Randeria (1995) :
[TABLE]
where =, =, is the -dependent current operator in the direction. Eq.(2) has been employed for metal-insulator transitions in the Hubbard model in many worksMa et al. (2018); Trivedi and Randeria (1995); Denteneer et al. (1999). A further way is to extract the spectral function by inverting the Laplace transform:
[TABLE]
in which can be achieved from the spatial Fourier transform of , and is solved by performing with a method of analytic continuation.
We are also concerned about the magnetic properties of the system by studying the AFM spin structure factorSorella et al. (2012); *PhysRevX.6.011029
[TABLE]
where represents the unit cell number of the lattice, is the component of the spin structure factor operator, and the angle brackets signify Monte Carlo simulations. To further study the nature of different stages of the system, we calculated the local moment byEuverte et al. (2013)
[TABLE]
III Results and discussion
We first calculated the temperature dependence of conductivity with increasing interaction for a fixed value =0.3. Basically, at low indicates an insulating phase, while at low corresponds to a metallic phase. In Fig. 2 (a), we can see that at =0.0 and =1.0, the curve shows a concave down tendency and approaches the origin as the temperature decreases. This kind of low behavior suggests that the system exhibits insulating behavior. Interestingly, the behavior becomes metallic as the on-site interaction increases to =1.8. A further increase in the interaction to =2.5 reinforces the metallic behavior, but the advance to =4.5 destroys the metallic state thoroughly and changes the system into a Mott insulator phase. When the staggered potential increases to =0.6 for =2.5, the increase in suppresses the metallic behavior and insulating state develops, which also implies that the value of =2.5 is not strong enough to drive the metallic phase for =0.6.
We present more data in Fig. 2 (b) to emphasize the contrast and highlight the effect of the staggered potential. For =0, at =1.0 and =3.5, increases as is lowered. When we calculate larger values (=4.0 and =4.5), decreases as is lowered, which shows insulator behavior when the staggered potential is absent. For the =4.0 case, the insulating phase at =0 is displaced by a metallic phase when =0.6.
Fig. 2 reveals the interesting phenomenon that the electronic correlation may drive a band insulator metallic, and at a larger interaction, there is a second transition from the metal to a Mott insulator. To further explore this issue, we plotted Fig. 3. Fig. 3(a) shows that the transition from metal to Mott insulator is restored at =3.9 for fixed =0.0. In Fig. 3(b), the curves intersect at two points, =1.4 and =4.2. In the range of 0<$$U$$<$$U_{c1} and U_{c2}$$<$$U$$<5.0, the conductivity values at higher temperature exceeds those of lower temperature for the same . The system maintains an insulating phase. The opposite situation emerges within the range of U_{c1}$$<$$U$$<$$U_{c2}. The conductivity increases as the temperature decreases, which demonstrates the metallic phase. The largest conductivity value for different values remains near =3.0 at =0.3. The two crosspoints represent the transitions from band insulator to metal to Mott insulator. Fig. 3(c) confirms these transitions with different =2.6 and =4.3 at =0.5, and the largest conductivity occurs at approximately =3.5.
Interestingly, however, the position of the largest conductivity moves to larger as increases, roughly following the law of =2.5+2, as that shown in Fig. 3(d). This result is quite different from that of the ionic Hubbard model on a square lattice, where the largest conductivity remains near =2, as one might expect from the =0 analysisParis et al. (2007). In the Hubbard model on a square lattice, the charge density wave and local moments are perfectly balanced on the special =2 line at =0, and hence, the system is most likely to be metallic. At =1, the AFM point also lies on that line, which is =0 at =0, but for a honeycomb lattice, the AFM point lies much higher above the line due to its Dirac fermion behavior at half filling. Therefore, perhaps the AFM point “pulls” the largest point up from =2, and this “pull”, also leaves a larger region of intermediate phase for interacting Dirac fermions with staggered potential.
To further support our analysis of the intermediate phases shown in Figs.1 to 3, we calculated the spectral functions for =0.5 and for three values of corresponding to the three phases. As shown in Fig.4, for =0.5 and =5.0, corresponding to the band insulating and Mott insulating phases respectively, shows a gap around =0, thus indicating an insulator. In contrast, for =3.5 is nonvanishing and shows a quasiparticle peak at =0, thus indicating a metallic phase. The results of spectral functions are consistent with the measurements of conductivity, as gets larger for =0.5, the spectral function shows an interaction-induced closing of a band gap and a subsequent opening of a Mott gap.
Fig. 5 provides the finite-size scaling results of the AFM structure factor on lattices of size =3, 6, 9, 12, 15. By extrapolating the data to the thermodynamic limit, we estimate that the critical point for is =4.04.3 when the band gap =0, which is consistent with the previous studies of AFM orderSorella et al. (2012); *PhysRevX.6.011029. As we can see, suppresses the AFM structure factor, while plays an opposite role. An increase in the on-site interaction to =4.5 AFM order for =0.3, and a further increase to =5.0 enables all calculated values to exhibit AFM order.
Local moment formation has been reported to happen discontinuously for the onset of Mott insulator behavior Euverte et al. (2013); Sentef et al. (2009); Tahara and Imada (2008). Fig.6 (a) shows on for a range of fixed values of for =12 and =12. We find that the Mott metal-insulator transition is associated with a maximum value in . The value of where the maximum in appears is very close to the characterized by the conductivity shown in Fig.1. Although the error bars make the value of transition point somewhat uncertain, the discontinuity of is an indicator of the phase transitions. The behavior of versus is given in Fig.6 (b). At small interaction (=2.0, 3.0) the system undergoes a very robust metallic phase, then immediately changes into a band insulator phase as gets larger, thus distinguishing the corresponding lines from the large ones, as at large the phase transition of the system falls near the Mott metal-insulator boundary.
The local moment is related with the double occupancy of =, and the evolution of may explain the phase transitions of the system to some extent. At ==0 and \Delta$$>0, the system is an insulator with some double occupancy because of the gap. If is increased, the charges would rearrange to reduce the double occupancy. If increased enough, eventually there is an elimination of double occupancy, one spin 1/2 occupies each site, and the model has long-range AFM order.
We also compared our results with those of an experimental studyMesser et al. (2015), in which their MI correlation image lies in our MI phase and CDW ordered phase lies in our band insulator (BI) phase. Their noise correlation image of the metal phase is a bit outside our metallic phase (which ends at =3.9). However, our value for , based on the conductivity, is consistent with the very precise simulations of Otsuka et al.Sorella et al. (2012); *PhysRevX.6.011029 which finds =3.9 based on the magnetic structure factor. Beyond this consistency, we obtain a phase diagram with numerically exact phase boundaries, making our results surpass the experimental work.
To make a more direct comparison with the experimentsMesser et al. (2015), the double occupancy as a function of the on-site interaction for various values of is shown in Fig. 7. For strong attractive interactions, a large fraction of doubly occupied sites is observed and it continuously decreases as increases. In the weak repulsive interactions where , the double occupancy continues to be large, as is that expected for the band insulator. For the strong interaction region of , the double occupancy tends to vanish for the largest , consistent with a Mott insulating phase. These trend and physical characteristics of double occupancy basically match the experimental results shown in Fig.2 of Ref. Messer et al. (2015).
IV Conclusions
Between the band insulating phase and Mott insulating phase, we find a metallic phase in the ionic Hubbard model on a honeycomb lattice. On a square latticeGarg et al. (2006); Paris et al. (2007); Bouadim et al. (2007), Coulomb interactions produce an AFM insulating phase at infinitesimal , and the competition with the staggered potential induces a metallic phaseBouadim et al. (2007). By contrast, on the honeycomb lattice, a Mott insulator transition occurs at finite even without staggered potential, and a rather wide region of and for the metallic phase is found.
In summary, we studied the ionic Hubbard model on a honeycomb lattice by a determinant quantum Monte Carlo method. We found that the intermediate phase between the two insulating phases was metallic. The staggered potentials drive the metallic phase to the band insulating phase, while the effect of the Coulomb repulsion was quite different. The Coulomb repulsion first drives the metallic phase into a nonmagnetic Mott insulating phase and then to the antiferromagnetic Mott insulating phase. As the Coulomb repulsion increases, the critical value of the staggered potential increases, suggesting competition between the two energy scales. However, along the Mott metal-insulator transition line, the effect of staggered potential is weak as the electrons are localized in the Mott insulating phase. Compared with previous theoretical proposals on some other models, our extensive numerical studies succeed in achieving a complete phase diagram, where the intermediate state is robust and occupies a large part of the phase diagram, and thus it should be more feasible to be detected experimentally.
Acknowledgements.
We thank Richard T. Scalettar and James E Gubernatis for many helpful discussions. We thank Hui Shao, Jinghua Sun and Tang Ho Kin for helping us to calculate spectral functions. This work was supported by NSFC (No. 11774033, No. 11974049, No. 11974048 and No. 11504023) and Beijing Natural Science Foundation (No. 1192011). We acknowledge the support of HSCC of Beijing Normal University, and some numerical simulations in this work were performed on Tianhe in Beijing Computational Science Research Center.
Appendix A Fintie size effect
To make the phase diagram more accurate and convincing, order parameters computed on finite-size lattices must be extrapolated to the thermodynamic limit. The finite-size effect on the spin structure factor was carefully examined in the manuscript. Here we discuss the size effect of the DC conductivity.
In Fig. 8, we plot the conductivity as a function of the temperature for lattices up to at metallic states (a) and insulating states (b). Both Fig. 8 (a) and (b)(=1.0) indicate that the lattice size has a distinct influence on the conductivity for . This result is predictable because the finite-size effects have remarkable impact on weak coupling. At and , there is an increase in with decreasing for the lattices that we have studied. Additionally, the metallic behavior weakens as the lattice size is increased. Although decreases with increasing lattice sizes, the signature of metallic behavior is unchanged. At and , the system shows an insulating behavior at low temperature, and results on larger lattice sizes reconfirm this behavior. At larger interaction strength as for the insulating states, the conductivity is almost independent of the lattice size.
These findings are consistent with the consensus that in a gapped system, one expects finite-size effects to be much smaller than in a metallic one. Because our focus is to discern the insulating phase, the data suggest that the lattice is large enough to be simulated for , and we could ascertain the insulating phase at low T.
Appendix B Zero temperature limit
The numerical method we employ, the finite temperature determinant quantum Monte Carlo (DQMC) method, can only calculate results at finite temperature. But we have ensured that the numerical results are converged at sufficiently low temperature, low enough to be regarded as the ground-state properties (zero temperature). Fig.9 is the detailed example. In Fig.9, we plot the AFM spin structure factor as a function of the inverse temperature under circumstances of different lattice sizes and staggered potentials.
The figure shows that the AFM order increases as the temperature is lowered. When drops below a lattice-dependent temperature, saturates and gradually level off, barely no dependence within acceptable statistical errors. So we reasonably conclude that the physical observable has reached the ground state if its value is convergent below some . In this work, we evaluated the limit to be for the spin structure factor. All the data presented in our paper were acquired at temperatures lower than . Therefore, the results we obtained can be regarded as the ground-state properties under the corresponding conditions. Considering this, we can do the size scaling.
Appendix C The DC conductivity formula
In our work, the low temperature behavior of DC conductivity is used to distinguish metallic or insulating phases. We implemented the approach proposed in the work by Trivedi et al.Trivedi et al. (1996), which is based on the following argument. From the fluctuation-dissipation theorem,
[TABLE]
where is the current-current correlation function. While could be computed by a numerical analytic continuation of data, we instead here assume that below some energy scale . Provided the temperature T is sufficiently lower than , the above equation simplifies to
[TABLE]
which is Eq. (2) in the manuscript.
It has been noted that this approach may not be valid for a Fermi liquidTrivedi et al. (1996). In this situation, the characteristic energy scale is set by . The requirement will never be satisfied. However, in the system we studied, the energy scale is set by the temperature-independent staggered potential strength so that Eq.(7) is valid at low temperatures.
It is worth mentioning that the authors of Refs.Trivedi et al. (1996); Ma et al. (2018) that used the DC conductivity formula were simulating disordered systems, whereas it still applies to a system without disorder. In Eq.(6), if we set to its largest value, , where is large, then the factor dies off rapidly and only small values contribute to the integral. More specifically, we expect a low frequency behavior where are important, and we can substitute this in and do the integration. This yields the approximate formula Eq.(7). Notice that for a Fermi liquid, so it is impossible to satisfy . Hence we cannot use Eq.(7). We are only safe if there is another scale (scattering mechanism), like disorder, which sets . For example, if we have disorder of strength and we can reduce to the point where and use Eq.(7).
Now, even though we do not have disorder, we do have an energy scale associated with the staggered potential . If we Fourier transform this, it does scatter the fermions but only between and , whereas a random potential mixes all .
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Novoselov et al. (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438 , 197 (2005) .
- 2Zhang et al. (2005) Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Nature 438 , 201 (2005) .
- 3Castro Neto et al. (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81 , 109 (2009) . · doi ↗
- 4Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82 , 3045 (2010) . · doi ↗
- 5Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83 , 1057 (2011) . · doi ↗
- 6Peres (2010) N. M. R. Peres, Rev. Mod. Phys. 82 , 2673 (2010) .
- 7Mourik et al. (2012) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336 , 1003 (2012) . · doi ↗
- 8Grover et al. (2014) T. Grover, D. N. Sheng, and A. Vishwanath, Science 344 , 280 (2014) . · doi ↗
